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1 . INTRODUCTION 

The  Captive  Trajectory  Yawmeter  System  (CTYS)  is  the  method  used  in  Aero- 
ballistics  Division,  WSRL  for  simulating  store  separation  from  aircraft. 

The  initial  plan  was  reported  in  reference  1,  preliminary  validations  are 
described  in  reference  2,  and  progress  reports  were  given  in  references  3 and  4. 
This  report  describes  the  mathematical  background  and  the  software  which  drives 
the  system. 

Section  2 provides  the  details  of  all  coordinate  systems  used,  defines  Euler 
angles  and  derives  the  equations  of  motion  for  the  six-degree-of-freedom  program. 
Section  3 outlines  the  methods  for  determining  the  forces  and  moments  on  the 
store,  while  Section  4 gives  details  of  the  system  and  software. 

The  CTYS  resides  on  a UNICHANNEL-15  (PDP-15/PDP-11)  dual  processor  computer 
which  is  attached  to  Aeroballistics  Division's  SI  wind  tunnel,  a continuous  flow 
tunnel  which  can  operate  in  a speed  range  from  about  100  m/s  to  Mach  1 and  from 
Mach  1.4  to  Mach  2.8.  A full  description  of  the  system  can  be  found  in 
reference  5 and  its  associated  references. 

In  view  of  the  small  size  (0.4  m x 0.4  m)  of  the  working  section,  the 
traditional  method  of  using  a model  of  the  store  mounted  on  a sting  has  been 
replaced  by  the  technique  of  measuring  the  flow  at  twelve  stations  along  the 
store  centre-line  for  each  trajectory  point  computed.  Thus,  the  store  loads 
are  derived  from  measurements  of  the  flow  field  and  are  not  directly  measured. 
Aircraft  models  are  typically  l/50th  scale. 

The  flow  measurements  are  made  by  a yawmeter  probe  mounted  on  a computer- 
controlled  four-degree-of-freedom  traverse  rig,  the  roll  axis  being  inoperative 
for  this  application.  The  probe  measures  pitot  pressure,  manifold  pressure  and 
two  differential  pressures,  which  are  sufficient  to  determine  the  flow  properties 
at  each  measurement  station. 

The  CTYS  then  determines  forces  and  moments  by: 

(a)  selecting  the  most  appropriate  result  from  data  which  have  been  previously 
acquired  on  an  approximately  l/8th  scale  model  of  the  store  in  uniform 
flow,  and 

(b)  calculating  increments  to  the  store- loading  due  to  the  known  non-uniform 
flow,  as  measured  on-line  by  the  probe. 

Section  3 gives  a full  account  of  this  method. 

In  other  respects  the  CTYS  is  analogous  to  the  Captive  Trajectory  System, 
being  simply  a point  prediction  technique  using  a fourth-order  Runge  Kutta 
integration  in  time.  One  significant  feature  of  the  method  of  load  estimation 
is  that,  as  the  store  moves  away  from  the  aircraft,  the  flow  non-uniformities 
decrease  and  the  corrections  to  the  measured  coefficients  tend  to  zero.  At  this 
stage  the  aircraft  does  not  influence  the  trajectory  any  more  and  the  trajectory 
can  then  be  continued  to  impact,  or  as  far  as  desired,  without  further  tunnel 
information. 

The  time  required  to  complete  an  average  trajectory  is  about  30  min.  Although 
the  system  has  been  streamlined  as  much  as  possible,  there  are  many  inherent 
delays  such  as  traverse  rig  travel  time  and  yawmeter  probe  settling  delay,  the 
latter  being  approximately  5 s.  Many  accesses  to  the  various  disk  data  stores 
must  be  made  for  each  trajectory  point,  but  these  are  performed  using  Direct 
Access  Input/Output,  mainly  while  the  rig  is  traversing  to  the  next  point,  and 
therefore  do  not  generally  delay  the  process. 

Results  to  date  are  promising  and  it  is  expected  that  full  scale  trajectories 
can  be  simulated  quite  successfully. 
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2.  MATHEMATICAL  ANALYSIS 

2.1  Definition  of  coordinate  systems  and  Euler  angles 

The  Captive  Trajectory  Yawmeter  System  uses  two  distinct  sets  of  coordinate 
axes.  The  first  set  refers  to  the  real  world  and  consists  of  earth  axes, 
aircraft  axes  and  bomb  axes  as  shown  in  figure  1. 

The  aircraft  is  assumed  to  move  in  a circle  in  the  X^Z^  plane  at  a constant 

tangential  velocity  and  with  centripetal  acceleration  pg.  Its  initial 

attitude  is  specified  by  a dive  angle  (DIVE)  and  an  angle  of  attack  (ATTACK). 
Also,  since  the  Aircraft  Z axis  (Z^)  points  in  a downward  direction  and  the 

Earth  Z axis  (Z^)  points  upwards,  the  roll  of  the  aircraft  with  respect  to 

earth  axes  is  it  radians  and  its  pitch  is  given  by 

APITCH  = -(DIVE-ATTACK). 

The  aircraft  maintains  this  trajectory  throughout,  where  only  the  dive 
angle  changes  with  time.  The  value  of  M in  the  centripetal  accelerations 
may  also  be  negative. 

The  CTYS  also  defines  a set  of  wind  tunnel  coordinates  as  shown  in 
figure  2.  The  flow  (VJ  is  parallel  to  the  tunnel  walls  whereas  the  air- 
craft and  bomb  axes  are  inverted,  the  aircraft  being  fixed  to  the  tunnel  wall 
at  an  angle  of  attack;  the  bomb  position  is  defined  by  a yawmeter  probe 
attached  to  a four-degree-of-freedom  traverse  rig.  The  wind  tunnel  Z axis 
(Z-j.)  points  downwards  thus  providing  a certain  parity  with  the  real  world 

system  of  figure  1. 

Throughout  the  analysis,  the  orientation  of  one  axis  system  with  respect 
to  another  is  described  by  Euler  angles  for  roll,  pitch  and  yaw.  For 
example,  to  convert  from  aircraft  axes  to  bomb  axes: 

(i)  Roll  through  an  angle  *f>  about  the  aircraft  X axis  (X^)  to  give  a 
new  system  (X'.Y'.Z'),  where  X'  = XA.  This  coordinate  trans- 
formation is  described  by  the  3x3  matrix 


C, 


10  0 
0 cos  p sin  'P 

0 -sin  $ cos 


(1) 


(ii)  Then  pitch  through  an  angle  0 about  the  Y'  axis  to  give  the  system 
(X",Y",Z"),  where  Y"  = Y'.  The  corresponding  transformation 
matrix  is 


1 

f cos  0 

0 

-sin  0 \ 

C2  = 

0 

1 

0 

(2) 

\ 

ysin  0 

0 

cos  0 J 

(iii)  Finally,  yaw  through 

an  angle  P 

about  the  Z" 

axis  to  give  bomb  axes 

w^ere  = Z".  This  corresponds  to  the  transformation 
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cos  P sin  P 0 

C3  = | -sin  P cos  P 0 

0 0 1 


The  total  transformation  is  thus  represented  by  the  direction  cosine  matrix 
C = C3C2C1 , which  is  given  by 


C = 


'cosO  cosp 
-cos 0 sin/3 
sin0 


sin^  sin#  cos P + cos#  sin P 
-sinp  sin#  sin/3  + cos<P  cos P 
-sin^  cos# 


-cos<£  sin#  cos P + siro/>  sin/3 \ 
cost?  sin # sin/3  + sin^  cos/3  I 
cos cosO  / 

(4) 


Thus,  if  a and  b are  vectors  in  aircraft  axes  and  bomb  axes  respectively, 
then 


b = Ca 

2.2  The  relation  between  Euler  rates  and  bomb  angular  rates 


(S3 


Let  p,q  and  r be  the  angular  rates  of  rotation  of  the  rotating  axis 
system  about  its  own  X,Y  and  Z axes  respectively.  To  calculate  the  relation- 
ship between  these  angular  rates  and  the  Euler  rates  , 9 and  P , it  is 
necessary  to  use  the  results  of  the  previous  section. 

Firstly,  a roll  of  <p  about  the  aircraft  X axis  gives  rise  to  a rate  vector 
(^,0,0)  in  aircraft  axes.  This  produces  a component 


co 


- C3  C2  C! 


(6) 


in  the  rotating  axis  system,  where  the  matrices  C, , C2  and  C3  are  defined  by 
equations  (1),  (2)  and  (3)  respectively. 

Secondly,  a pitch  of  6 about  the  Y'  axis  (see  Section  2.1)  produces  a 
rate  component 


co  2 - C3  C2 


(7) 


in  the  rotating  axis  system. 

Finally  a yaw  of  P about  the  Z"  axis  produces  a rate  component 


co 


= C3 


(8) 


in  the  rotating  axis  system. 

Summing  the  contributions  of  equations  (6),  (7)  and  (8),  we  obtain 
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p = 6 sin  0 + ' fi  cos  0 cos  0 

q = 0 cos  0 - ¥>  cos  0 sin  0 

and 

r = 0 + $ sin  0 . 

Conversely, 

'f  = (p  cos  0 - q sin  0}/cos  0 
0 = p sin  0 + q cos  0 

and 


(9(a)) 

(9(b)) 

(9(c)) 

(10(a)) 

(10(b)) 


0 = r - tan  0 (p  cos  0 - q sin  0).  (10(c)) 

Note  that  equations  (10)  have  a singularity  at  0 = j.  This  fact  will  be  of 
significance  later. 

2.3  Rate  of  change  of  the  Euler  angle  transformation  matrix 

It  will  be  useful  to  calculate  the  rate  of  change  of  the  direction  cosine 
matrix  C depicted  by  equation  (4) . 

Let  F be  a fixed  axis  system  and  B be  an  axis  system  which  is  rotating  at 
a rate  to  with  respect  to  F. 

Let  b be  a vector  in  B whose  coordinates  in  the  fixed  system  F are  given 
by  f.  Then 


b = Cf  (11) 

where  C is  a 3 x 3 orthogonal  direction-cosine  matrix. 

Now,  it  is  well  known  that 


£ = b + to  x b. 


(12) 


the  second  term  deriving  from  the  rotation  of  axis  system  B.  Now  since  to 
is  a vector  in  F,  then  we  can  define  a vector  to  ' = (p,q,r),  where  p,  q and 
r are  the  angular  rates  of  rotation  of  B about- its  own  X,  Y and  Z axes 
respectively.  Further,  since  to  ' is  a vector  in  B,  then 


to  ' = Cto  . 


Equations  (12)  and  (13)  then  give 


f = b + C'to'xb. 


(13) 


l 


(14) 


n 


Now  w ' x b can  be  expressed  as  the  matrix 

/ 0 r 

fi  = j -r  0 

\ q -p 

Hence,  equations  (14)  and  (15)  yield  the  relations 

f = b - C"1  flb.  (16) 

To  calculate  C,  the  rate  of  change  of  the  direction  cosine  matrix  C,  let  b 
be  a vector  which  is  fixed  in  the  rotating  system  B,  so  that 

b = 0 (17) 

and 

f = -C"1  n b . (18) 
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product  b,  where 


(15) 


Differentiation  of  equation  (11)  with  respect  to  time  gives 

b = Cf  + Cf . (19) 

Then  substitution  of  equations  (11),  (17)  and  (18)  into  equation  (19)  yields 


That  is. 


0 = Cf  + C(-C~‘  n b) 

= Cf  - ft  C f. 


(C  - ftc)f  = 0 


(20) 


Since  f is  an  arbitrary  vector,  we  can  deduce  that 

c = n C. 


(21) 


2.4  Trajectory  of  aircraft 

As  noted  in  Section  2.1,  the  aircraft  describes  a circular  trajectory 
which  is  either  convex  to  the  origin  (p  positive)  or  concave  to  the  origin 
(p  negative).  The  centripetal  acceleration  (pg)  is  assumed  to  include  the 
gravitational  component. 

Thus,  at  all  times 


= Mg 


(22) 
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- 6 - 


" dt  (DIVE) ' 


(23) 


where  p is  the  radius  of  curvature  and  DIVE  is  the  dive  angle,  a function 
of  time. 

Eliminating  p from  equations  (22)  and  (23)  gives 

DIVE (t)  = DIVE  - Pgt/V  . (24) 

o 00 

Further  analysis  then  yields  the  coordinates  of  the  aircraft  in  earth  axes 


as  follows 

x ' x0  * 7F  “S(DIVE0  - lfr)sin  ($;) 

y = 0 (25(b)) 

2 ' zo-irg  sin(DIVE°  ‘ “)sinfe)'  (25(c)) 

When  p is  small,  equations  (25)  reduce  to 

x = xq  + V^t  cos(DIVE^)  (26(a)) 

Y = 0 (26(b)) 

z = Zo  - ^ sin(DIVEo).  (26(c)) 


The  case  p = 0,  of  course,  corresponds  to  a straight  line  trajectory. 

From  figure  1,  it  is  clear  that  the  components  of  aircraft  acceleration 
in  aircraft  axes  are  constant  and  are  given  by 

xA  = pg  sin (ATTACK)  (27(a)) 


yA  = 0 (27(b)) 

zA  = -Pg  cos (ATTACK).  (27(c)) 


Also,  the  Euler  angles  describing  the  orientation  of  the  aircraft  with 
respect  to  earth  axes  are 

AROLL  = n (28(a)) 


APITCH 


-(DIVE-ATTACK) 


128(b)) 


Substitution  of  these  angles  into  the  direction  cosine  matrix  and 
application  of  the  transformation  x^  = 1 x^  gives  the  components  of 

aircraft  acceleration  in  earth  axes,  namely. 


xE  = Pg  sin(DIVE)  (29(a)) 

XE  = 0 (29(b)) 

zE  = (*g  cos  (DIVE).  (29(c)) 

2.5  Derivation  of  the  equations  of  motion 

We  must  now  derive  the  six  equations  of  motion  which  describe  the  motion 
of  the  bomb  relative  to  the  aircraft. 

Let 

(i)  X be  the  position  vector  of  the  bomb  and  X*  be  the  position  vector 
of  the  aircraft  with  respect  to  some  axis  system, 

(ii)  Subscripts  E and  A refer  to  earth  and  aircraft  axes,  and 

(iii)  CA  be  the  3x3  direction-cosine  matrix  describing  the  orientation 

of  the  aircraft  relative  to  earth  axes. 

Then 

(i)  ^ is  the  position  vector  of  the  bomb  relative  to  earth  axes, 

(ii)  X£  is  the  position  vector  of  the  aircraft  relative  to  earth  axes, 
and 

(iii)  is  the  position  vector  of  the  bomb  relative  to  aircraft  axes. 
Thus, 


(30) 


(31) 


(32) 
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where  u would  describe  the  velocity  of  the  bomb  centre  of  gravity  with 
respect  to  aircraft  axes  if  were  constant,  that  is,  if  the  aircraft  were 

not  rotating  about  its  own  centre  of  gravity. 

Then  equations  (33)  and  (30)  in  equation  (32)  give 

h = H + n A Xa  (34) 

and  ^ of  equation  (33)  yields 

“ - CA  X£  - CA  X*  + nA  u.  (35) 

Now  Xg  is  the  acceleration  of  the  bomb  with  respect  to  earth  axes,  expressed 

in  earth  axes.  Thus,  since  earth  axes  represent  an  inertial  system,  we  can 
write 

Xe  = £5/1",  (36) 

where  F^  are  the  external  forces  on  the  bomb,  expressed  in  earth  axes,  and  m 
is  the  mass  of  the  bomb. 

Then,  premultiplication  of  equation  (36)  by  the  3x3  matrix  C expresses 

t\ 

the  forces  in  aircraft  axes;  thus 

CA  h = CA  V"1  = Vm‘ 

Also,  since  X£  represents  the  acceleration  of  the  aircraft  relative  to  earth 
axes,  then  CA  represents  the  same  acceleration  but  now  expressed  in 
aircraft  axes.  Thus,  by  equations  (27),  we  have 


C.  X*  = 
A —E 


pg  sin (ATTACK) 


-pg  cos (ATTACK) 


Finally,  substitution  of  equations  (37)  and  (38)  into  equation  (35)  gives 


u = FA/m  - 


pg  sin (ATTACK) 


,-pg  cos  (ATTACK) 


+ SI  u. 
A — 


Equations  (34)  and  (39)  are  thus  the  two  first-order  differential  equations 
which  can  be  integrated  to  give  the  motion  of  the  bomb  centre  of  gravity 
relative  to  aircraft  axes. 

To  derive  the  equations  of  motion  describing  the  attitude  of  the  bomb 
relative  to  aircraft  axes,  we  define  li  the  angular  momentum  vector  as 

“O 
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^B 


(40) 


where  I I^y  and  are  the  moments  of  inertia  of  the  bomb  about  its  own 

principal  axes  and  p,  q and  r are  the  angular  rates  of  the  bomb  about  these 
same  X,  Y and  Z axes  respectively. 

The  inertia  tensor  is  not  a constant  if  it  is  calculated  in  the  space 
coordinate  system  (earth  axes)  because,  in  this  system,  the  values  of  the 
coordinates  change  with  time.  In  the  moving  (bomb)  system,  however,  the 
inertia  tensor  is  constant.  For  this  reason  it  is  advantageous  to  work  in 
the  bomb  system. 

If  h^  is  the  vector  representing  the  angular  momentum  of  the  bomb  in  earth 
axes,  then 


^B 


= CB  Mg* 


(41) 


where  CD  is  the  direction-cosine  matrix  which  describes  the  orientation  of  the 

D 

bomb  relative  to  earth  axes, 
of  equation  (41)  gives 

(42) 

(43) 

Now,  since  earth  axes  represent  an  inertial  system,  we  can  write 


hfi  = CB  + CB  Mg- 


Then  equations  (21),  (41)  and  (42)  give 


h 


CB^  + nB  V 


(44) 


where  is  the  vector  representing  external  moments  (torque)  on  the  bomb, 
expressed  in  earth  axes. 

Thus,  premultiplying  equation  (44)  by  CD  expresses  the  moments  in  bomb 

D 

axes,  namely 


C 


B 


Mb- 


(45) 


Finally,  equations  (43)  and  (45)  yield 


hfi  Mb  + ^ B 


(46) 


Equations  (40)  and  (46)  thus  enable  the  calculation  of  the  bomb  angular 
rates  p,  q and.r..  In  most  cases,  equations  (10)  would  then  give  the  Euler 
angular  rates  d and  $ which  could  then  be  integrated  to  give  the  Euler 
angles  0 and  (1  as  functions  of  time. 
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However,  equations  (10)  have  a singularity  at  0 = 7r/2  and  thus,  for 
values  of  0 in  the  neighbourhood  of  7r/2,  truncation  errors  will  occur. 

It  is  therefore  safer  to  use,  not  the  Euler  angles  themselves,  but  the 
direction-cosines  of  these  angles  as  dependent  variables. 

Thus,  let  C be  the  direction-cosine  matrix  which  represents  the  orient- 
ation of  the  bomb  relative  to  the  aircraft.  If  b is  a vector  fixed  in 

bomb  axes  and  a is  a vector  fixed  in  aircraft  axeJ,  then 

b = Ca.  (47) 

Further,  if  e is  a vector  in  the  earth  axis  system,  we  have 


n = ci£ 


(48) 


and 


a = CAe. 


(49) 


Simple  manipulation  then  gives 


CCA  * CB 


(50) 


or 


c = cBcA\ 


(51) 


of  equation  (50)  yields 


CC.  + CCS 
A A 


= C, 


(52) 


Equation  (21)  applied  to  equation  (52)  thus  gives 


CCA  + CfiACA  = ^BCB- 


(53) 


Post-multiplying  equation  (53)  by  C^'1  and  substituting  equation  (50)  finally 
gives 


c = ft  n c - c n 

B A 


(54) 


The  3x3  matrix-equation  (54)  thus  represents  9 first-order  differential 
equations  for  the  calculation  of  the  elements  of  the  orientation  matrix  C. 
The  procedure  is  therefore  to 

(i)  Integrate  equation  (46)  to  calculate  h_. 

— B 

(ii)  Use  equation  (40)  to  give  p,  q and  r. 

(iii)  Substitute  in  equation  (15)  to  give  H 
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(iv)  Calculate  H A from  the  known  aircraft  angular  rates  p^,  q^  and  r^ 

(v)  Substitute  into  equation  (54)  and  integrate  it  to  yield  the 
3x3  matrix  C. 

If  so  desired,  the  Euler  angles  associated  with  C can  be  calculated  by 
the  following  formulae,  which  are  derived  from  equation  (4):- 


0 = tan"1 (C31/VC32  + Ci3) 


sin  P = -C2 1 /cos  9 


cos  P = C, , /cos  9 


(55(a)) 


(55(b)) 


(55(c)) 


sin  <p  = -C.12  /cos  9 


(55(d)) 


cos  <p  = C33/COS  9 


(55(e)) 


Equation  (55(a))  assumes  that  9 is  acute. 

However,  as  noted  in  the  next  section,  it  is  usually  more  meaningful  to 
use  projected  angles  rather  than  Euler  angles. 

2.6  Transformations  between  Euler  angles  and  projected  angles 

Euler  angles  are  always  difficult  for  the  experimenter  to  visualize. 

For  instance,  in  defining  the  orientation  of  the  bomb  relative  to  aircraft 
axes,  it  is  extremely  difficult  to  relate  the  bomb's  attitude  to  the  given 
Euler  angles.  It  is  therefore  convenient  to  express  the  attitude  in  terms 
of  projected  angles. 

Thus,  let  •fi* , 9*  and  P*  be  the  projected  roll,  pitch  and  yaw  of  the  bomb 
relative  to  aircraft  axes,  that  is,  the  roll,  pitch  and  yaw  angles  of  the 
image  of  the  bomb  when  projected  upon  the  aircraft  YZ,  XZ  and  XY  planes 
respectively. 

To  find  projected  pitch  ( 9* ) and  projected  yaw  (/?*)  it  is  necessary  to 
find  the  coordinates  of  the  unit  vector  in  the  direction  of  the  bomb's  nose 
(X  axis)  as  represented  in  aircraft  axes.  This  is  given  by 


where  C is  the  matrix  describing  the  orientation  of  the  bomb  relative  to 
aircraft  axes. 

Since  C is  an  orthogonal  matrix,  then  equation  (56)  gives 


XM  = C. 


(57(a)) 


YN  = C, 


(57(b)) 


ZN  = C, 


(57(c)) 


9*  = tan"1 (-ZN/XN)  = -tan"1 (C, 3/C, , ) 
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and 


P*  = tan'1 (YN/XN)  = tan"' (C, 2 /C, , ) . 


(59) 


Similarly  the  projected  roll  (>p*)  can  be  found  from  the  coordinates  of 
the  unit  vector  in  the  direction  of  the  bomb's  reference  fin  (Z  axis)  as 
represented  in  aircraft  axes,  given  by 


= C"M  0 


/C3A 

C32 

\c3J 


sin  0 
-sin  p cos  9 
cos  p cos  6/ 


Since  roll  is  always  the  first  of  the  Euler  angles,  projected  roll  (p*)  and 
Euler  roll  (tp)  are  equal,  giving 


sin  p*  = -C32/COS  0 


(60(a)) 


and 


cos  p*  = C33/COS  d (60(b)) 

Differentiation  of  equations  (60)  with  respect  to  time  shows  that 


'P*  — p = (C33  C32  - C32  C3  3 ) / (C3  2 + C33) 


(61) 


Projected  pitch  and  yaw  rates  are  obtained  by  differentiating  equations  (58) 
and  (59)  to  give 


and 


- (Ci  1 Ci  3 - Ci  3 C,  1 )/(C? , + Ci  3) 


P*  = (C,2  C, , - C, , C,  2 )/ (C? . + Ci  2 ) • 


(62) 


(63) 


Further,  equation  (4)  can  be  used  to  substitute  for  the  elements  of  C to 
give  the  projected  angles  in  terms  of  the  Euler  angles,  namely:- 


and 


if>*  = ^ 


tan  0*  = tan  0 cos  f>  - sin  ^ tan  P/cos  0 


tan  P*  = tan  0 sin  + cos  f>  tan  P/ cos  9 


(64(a)) 

(64(b)) 

(64(c)) 


The  inverses  of  equations  (64)  are  then 


ip  = 


(65(a)) 
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tan  0 = tan  0*  cos  <p*  + tan  P*  sin  <p*  (65(b)) 

and 

tan  P = cos  0 (tan  P*  cos  y?*  - tan  0*  sin  y>*).  (65(c)) 

2.7  Relation  between  Projected  Angular  Rates  and  Bomb  Angular  Rates 

In  Section  2.2,  p,  q and  r were  defined  as  the  angular  rates  of  rotation, 
relative  to  fixed  axes,  of  the  rotating  system  about  its  own  X,  Y and  Z axes 
respectively.  Thus  we  here  define  p,  q and  r to  be  the  angular  rates  of 
the  bomb,  relative  to  earth  axes,  about  the  bomb's  X,  Y and  Z axes  respect- 
ively. 

If  p^,  qA  and  rA  are  the  angular  rates  of  rotation  of  the  aircraft, 

relative  to  earth  axes,  about  the  aircraft's  X,  Y and  Z axes  respectively, 
then  we  can  write 

£ = CEa  + d£,  (66) 

where 


C is  the  direction  cosine  matrix  describing  the  orientation  of  the  bomb 
relative  to  the  aircraft  and  d£  is  the  effect  of  the  rates  of  change  of  the 
Euler  angles  of  C. 

Now,  differentiating  equations  (65)  with  respect  to  time  gives 

f = (67(a)) 

sec2  0 0 = sec2  0*  cos  0*  + sec2  P*  sin  P*  + tan  P *p*/cos  0 

(67(b)) 

and 

sec2j 3 P = cosd(sec2j3*  costp*  P*  - sec20*  sin^*  0*  - tan 0 ^*)  - tan#  tan  P 0 

(67(c)) 


Then,  equations  (4),  (9)  and  (67)  give 
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dp  = 

C2 1 (C3  2 

sec2 

0*  0* 

- C33 

sec2  0*0*)  + cos  0 <p*/ cos  0 

(68(a)) 

dq  = 

Ci  1 (C33 

sec2 

0*0* 

- C32 

sec2  0*  0*) 

(68(b)) 

dr  = 

C,  1 (C2  2 

sec2 

0*  0* 

- C23 

sec 20*  0*) . 

(68(c)) 

Equations  (66)  and  (68)  therefore  provide  a means  of  calculating  bomb 
angular  rates  from  projected  angular  rates,  the  latter  being  specified  by 
the  experimenter  prior  to  the  running  of  the  Captive  Trajectory  Yawmeter 
System. 


3.  DETERMINATION  OF  THE  FORCES  AND  MOMENTS 

The  Captive  Trajectory  Yawmeter  System  determines  the  aerodynamic  forces  and 
moments  on  the  bomb  by 

(a)  selecting  the  most  appropriate  result  from  data  which  have  been  previously 
acquired  on  an  approximately  l/8th  scale  model  of  the  store  in  uniform 
flow,  and 

(b)  calculating  increments  to  the  store  loading  due  to  the  known  non-uniform 
flow,  as  measured  on-line  by  a yawmeter  probe. 

3.1  Store  characteristics  in  uniform  flow 

Tabulated  data  representing  the  characteristics  of  the  store  in  a uniform 
flow  are  required  for  this  procedure.  The  tables  reside  on  disk  and  consist 
of  ten  coefficients  as  a function  of  reference  Mach  number,  total  pitch  and 
total  roll. 

The  yawmeter  probe  is  programmed  to  take  readings  at  twelve  stations  along 
the  calculated  position  of  the  store  centre-line,  including  three  stations 
located  at  the  leading  edge,  mid-point  and  trailing  edge  of  the  mean  aero- 
dynamic chord  of  the  fins. 

A reference  point,  usually  the  station  at  the  mid-point  of  the  fins,  is 
selected  and,  from  these  reference  flow  conditions,  a Mach  number,  total 
pitch  and  total  roll  are  calculated.  These  quantities  are  then  used  to 
obtain  the  ten  uniform-flow  coefficients  from  the  disk. 

In  calculating  the  incidence  and  sideslip  upon  which  total  pitch  and  total 
roll  are  based,  some  allowance  must  be  made  for  the  curvature  of  the  stream 
past  the  fins.  Here  use  is  made  of  thin  aerofoil  theory  from  which  it  can 
be  shown  that  "the  lift  of  an  aerofoil  of  camber  7 and  incidence  a is  equal 
to  the  lift  of  a straight  aerofoil  at  incidence  (a  + 27)".  This  result  is 
easily  proved  using  page  91  of  reference  6. 

Consider  the  velocities  at  the  three  fin  stations  in  the  XZ  plane  of 
figure  3. 

The  slope  of  the  streamline  at  any  station  is  given  by 


^ = z/x  = W/U,  (69) 

where  U and  W are  the  measured  flow  velocity  components  in  the  X and  Z 
directions  respectively.  A parabola  is  fitted  through  the  streamline 
slopes  at  the  three  fin  stations  and  the  equation  of  the  streamline  is  thus 
obtained  by  integration  of  equation  (69)  with  respect  to  x. 
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The  equation  of  the  chord  (figure  3)  is  given  by 

zOO  = (zc/xc)  x-  (70) 

Thus  the  mean  height  is  given  by 

h = A/cos  a = (zfi  - zc  xB/xc)/cos  a (71) 

The  mean  chord,  quite  clearly,  is  given  by 

K = xc/cos  a (72) 

Thus,  the  camber  is  given  by 

y = h/K  = (zB  - zc  xB/xc)/xc.  (73) 

In  addition,  the  mean  incidence  is  given  by 

a = tan"' (zc/xc).  (74) 

Equations  (73)  and  (74)  thus  provide  the  means  for  calculating  the 
effective  incidence  ag  = (a  + 27). 

The  above  procedure  is  repeated  in  the  XY  plane  to  determine  the 
effective  sideslip  PQ  = (P  + 27). 

In  deriving  the  most  appropriate  uniform  flow  coefficients  from  the  disk 
data  store,  it  is  assumed  that  the  tail  loading  contained  in  these 
coefficients  is  representative  of  the  tail  loading  in  the  curved  flow  field. 

The  data  store  coefficients  are,  of  course,  based  on  free  stream  dynamic 
pressure.  Thus,  to  obtain  the  coefficients  appropriate  to  the  reference 
point,  the  data  store  coefficients  are  multiplied  by  reference  point  dynamic 
pressure  and  divided  by  free  stream  dynamic  pressure. 

3.2  Corrections  for  non-uniform  flow 

Since  the  interpolated  coefficients  of  the  previous  section  are  for  a 
uniform  flow,  the  flow  conditions  implied  at  each  store  measurement  station 
are  equal  to  those  at  the  store  reference  point.  However,  corrections  must 
be  applied  to  allow  for  changes  in  body  loading  forward  of  the  fins,  due  to 
the  differences  between  the  measured  flow  conditions  at  the  store  stations 
and  the  uniform  flow  implied  by  the  interpolated  coefficients.  * 

Allen  and  Perkins (ref. 7)  derive  an  expression  similar  to  that  derived  by 
Munk  for  the  potential  cross  force  on  slender  bodies,  given  by 


f = q ig  sin  2a,  (75) 

where  q is  the  dynamic  pressure,  and  a is  body-axis  incidence  relative  to 
the  free-stream  flow  direction. 

They  also  note  that,  from  the  work  of  Ward (ref. 8),  it  can  be  shown  that 
the  potential  cross-force  is  directed  midway  between  the  normal  to  the  axis 
of  revolution  and  the  normal  to  the  wind  direction. 
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This  applied  to  equation  (75)  gives  a cross-force  of 


ds 

ZP  = ^ dx  Sin  2a  COS  55  a 


(76(a)) 


YP  = ^ dx  Sin  ^ COS  * ^ 


(76(b)) 


for  the  Z and  Y potential  flow  contributions  respectively. 

A second  term  is  due  to  viscosity  and  can  be  written  as 

fv  = * CD  D q sin2  9t>  (77) 

c 

where  Cp  is  the  cross-flow  drag  coefficient, 

D is  the  cross-sectional  diameter  of  the  bomb  (i.e.  S = ttD2/4), 

i?  is  the  cross-flow  drag  proportionality  factor, 

0^  is  the  total  pitch  of  the  bomb  (see  figure  4), 

and  q is  dynamic  pressure  (=h  P I V|  2 ) . 

r)  arises  from  the  fact  that  the  ratio  of  the  cross-flow  drag  for  a finite 
length  cylinder  to  that  for  an  infinite  length  cylinder  is  less  than  one 
(see  reference  9) . 

Now,  since  f acts  in  the  total  Z axis  negative  direction,  the  components 

of  the  viscous  force  in  the  Z and  Y directions  will  be  given  by  (see 
figure  4) 


Zv  = -h  P V CD  D I V I 2 sin2  0T  cos 


(78(a)) 


Y v = -h  P T)  Cp  D I VI  2 sin2  0 T sin  ^ . 


(78(b)) 


But,  as  shown  in  figure  4, 


-I  V|  cos  0T, 


(79(a)) 


-I  Vl  sin  0T  sin 


(79(b)) 


W = -I  Vl  sin  0T  cos  1 Ap. 


(79(c)) 
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Thus,  equations  (79)  in  equations  (78)  yield 


ZV  = ^IC^DW  s/j2  * W2 


(80(a)) 


and 


Yv  = h P V cD  D V v^2  + W2 


(80(b)) 


Finally,  as  noted  in  reference  10,  allowance  must  be  made  for  buoyancy. 
This  results  in  the  components 


ZB  ' -V»'>SS 


(81(a)) 


and 


■v-'>  sf. 


(81(b)) 


where  is  free-stream  velocity  and  p is  local  air  density.  The  buoyancy 

term  arises  because  the  streamlines  at  any  point  are  curved.  This  is 

confirmed  by  measurements  taken  by  the  yawmeter  probe  along  the  calculated 
position  of  the  store  centre-line.  However,  buoyancy  is  not  expected  to 
contribute  a large  amount  to  the  overall  cross-force. 

The  Captive  Trajectory  Yawmeter  System  thus  defines  total  cross- 
forces of 


Z = Zp  + Zv  + ZB  (82(a)) 

and 

Y = Yp  + Yy  + Yb.  (82(b)) 

Uniform-flow  theory  says  that  flow  conditions  at  all  body  stations  are 
equal  to  those  at  the  reference  point.  Thus,  the  calculations  implied  by 
equations  (76)  to  (82)  are  repeated,  with  only  the  body  geometry  being 
allowed  to  vary  between  stations,  to  give  the  local  cross-force  predicted  at 
each  station  by  uniform  flow  theory,  namely  Z . 

Then  the  increment  in  cross-force  due  to  the  non-uniformity  of  the  flow 
at  each  station  is  given  by 


AZ  = Z - Z . 

u 


Similarly,  a pitching  moment  at  each  station  can  be  calculated,  with  resp-' 
to  an  origin  at  the  centre  of  gravity  of  the  store,  as 


i 

i 


M = Z(x  - xCG).  1 
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Thus  the  increment  due  to  non-uniformity  of  the  flow  is  given  by 

AM  = (z  - zu)  (X  - xCG). 

Finally,  this  analysis  can  be  applied  to  the  yaw  plane  to  yield  values  for 
Ay  and  yawing  moment  An. 

In  the  limit,  as  the  flow  becomes  uniform,  these  corrections  all  tend  to 
zero. 


4.  SYSTEM  DESCRIPTION 

The  software  for  the  Captive  Trajectory  Yawmeter  System  consists  of  an  overlay 
system  comprising  a resident  program  plus  several  links  which  overlay  each  other 
in  core.  This  has  been  necessary  because  of  the  limited  core  size  (24K  18-bit 
words)  of  the  UNICHANNEL-15  computer. 

All  relevant  variables  have  been  allocated  to  Blank  COMMON,  and  this  resident 
COMMON  area  is  written  on  disk  after  the  calculation  of  each  trajectory  point. 
This  feature  allows  the  trajectory  calculations  to  be  resumed  from  a previous 
trajectory  and  thus  provides  a useful  restart  facility. 

At  the  beginning  of  a new  trajectory,  the  COMMON  area  is  initiated  by  reading- 
in  a Constants  Block  from  the  disk.  The  Constants  Block  is  set-up  by  means  of 
a special  Editor  prior  to  the  running  of  an  experiment  and  resides  unchanged  on 
disk  until  it  is  required  to  change  any  of  the  parameters.  This  removes  the 
need  to  type-in  a large  amount  of  information  at  the  beginning  of  each  run. 

The  Constants  Block  includes  the  initial  position,  attitude  and  rates  of  the 
store  because  the  present  CTYS  begins  its  calculations  after  the  ejector  has 
operated  and  thus  the  effects  of  the  ejector  must  already  be  present  in  the 
Constants  Block  when  the  CTYS  begins  its  operation.  However,  the  ejector  phase 
of  the  CTYS  has  been  developed  separately (ref . 4)  and  work  is  currently  underway 
to  incorporate  it  as  an  integral  part  of  the  CTYS. 

4.1  Operation  of  the  Constants  Editor 

The  Constants  Editor  (EDCON)  either  sets  up  a new  file  or  allows 
modifications  to  be  carried  out  on  an  existing  file.  It  operates  in  an 
interactive  manner  by  prompting  the  operator  and  acquiring  answers  to  the 
questions  it  asks. 

Listings  of  EDCON  and  its  utility  routines  DECDE  and  IDECDE  are  given  in 
Appendix  I.  These  listings  specify  the  elements  of  the  arrays  T and  IT, 
which  together  comprise  the  COMMON  block  for  the  CTYS. 

The  loading  sequence  for  EDCON  is: 


$A  RKA  -5 
$GLOAD 

LOADER  V3A000 
> — EDCON  < ALTMODE  > 


When  EDCON  has  loaded,  it  announces  its  presence  and  asks  the  operator  to 
specify  a file  name.  This  file  name  will  denote  either  a new  file  to  be 
created  or  an  existing  file  to  be  modified. 

If  the  file  is  new,  EDCON  enters  an  Input  mode  and  acquires  the  values  of 
all  relevant  constants  by  means  of  a series  of  questions  and  answers. 

Each  prompt  by  EDCON  includes  the  units  in  which  the  answer  is  expected,  thus 
providing  a convenient  way  to  permanently  establish  the  units  of  each  input. 
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If  the  file  is  not  new,  EDCON  enters  an  Edit  mode  which  first  checks  for 
the  existence  of  the  file  and,  if  found,  opens  it.  It  then  prompts  the 
operator  with 


ALTER  PARAMETERS? 


A negative  answer  (N)  at  this  stage  would  simply  produce  a listing  of  the 
unmodified  Constants  Block.  However,  an  affirmative  (Y)  would  cause  a 
second  prompt  asking  the  operator  to  enter  N,  where  N denotes  the  following: 

(1)  Bomb  parameters 

(2)  Aircraft  parameters 

(3)  Probe  parameters 

(4)  Miscellaneous  information 

(5)  Measurement  Stations 

(6)  Diameters  at  Measurement  Stations. 

The  value  for  N directs  EDCON  to  the  appropriate  area  in  the  Constants 
Block,  where  it  begins  to  list  each  parameter  (just  as  it  did  in  Input  mode) 
together  with  the  current  value  of  that  parameter.  If  the  operator  wishes 
to  change  a parameter,  he  types  in  the  new  value.  A CARRIAGE  RETURN  will 
accept  the  old  value. 

The  operator  may  allow  listing  to  continue  to  the  end  of  the  section, 
changing  the  values  of  any  desired  parameters,  or  he  may  terminate  this 
phase  by  typing  CONTROL  P.  In  either  event,  EDCON  will  then  again  type 


ALTER  PARAMETERS? 


The  process  is  repeated  until  the  operator  finally  answers  N(0)  to  this 
question.  Then  EDCON 

(i)  calculates  the  aircraft  Euler  angles  and  their  rates,  relative  to 
earth  axes. 


(ii)  calculates  the  Euler  angles  defining  the  orientation  of  the  bomb 
relative  to  aircraft  axes. 


(iii)  determines  the  bomb's  maximum  cross-sectional  area,  and  its 
corresponding  diameter, 

(iv)  initializes  the  direction-cosine  matrix  defining  the  orientation 
of  the  bomb  relative  to  aircraft  axes,  according  to  equation  (4), 

(v)  calculates  the  angular  rates  of  the  bomb  as  described  in 
Section  2.7,  then 

(vi)  calculates  the  initial  values  of  bomb  angular  momentum  as  defined 
by  equation  (40), 

(vii)  finds  the  linear  accelerations  of  the  aircraft  relative  to  aircraft 
axes,  as  defined  by  equations  (27),  and  finally 

(viii)  calculates  the  free-stream  values  for  velocity,  velocity  of  sound 
and  ambient  temperature,  from  free-stream  Mach  number  and  the 
height  of  the  aircraft. 


EDCON  then  writes  the  file  on  disk  and  prints  the  values  of  all  entered 
constants  on  the  line  printer. 

The  process  of  entering  the  constants  is  straightforward  and  should  not 
cause  any  difficulty.  However,  it  should  be  remembered  that  the  roll, 
pitch  and  yaw  of  the  bomb,  together  with  their  rates,  are  entered  as 
projected  angles,  not  Euler  angles. 
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4.2  Operation  of  the  Main  Program 

The  Main  Program  (MAIN),  like  the  COMMON  area,  is  resident  in  core 
throughout  the  entire  trajectory. 

Upon  loading,  its  first  task  is  to  initiate  the  traverse  rig  console 
button.  This  primes  the  program  so  that  it  will  change  to  a "bomb-alone" 
trajectory  when  the  button  is  pressed,  namely,  when  the  operator  judges  that 
the  bomb  is  sufficiently  clear  of  the  aircraft  to  be  in  a nominally  uniform 
flow.  The  probe  is  then  no  longer  required  and  a trajectory  to  impact  is 
calculated. 

The  Main  Program  next  initiates  the  system  for  Direct-Access  Input/Output, 
which  permits  the  user  to  reference  directly  any  record  in  a file  without 
indexing  from  the  file's  beginning  up  to  the  desired  record.  This  feature 
is,  of  course,  essential  to  an  operation  such  as  the  CTYS,  where  very  large 
numbers  of  records  must  be  accessed  in  a random  manner. 

After  announcing  its  presence,  MAIN  prompts  the  operator  to  determine 
whether  the  trajectory  is  new  or  is  simply  the  continuation  of  a previous  one. 
If  the  trajectory  is  not  new,  the  system  reads  into  core  the  contents  of 
file  'TFILE  CTS',  the  last  image  of  the  COMMON  area  calculated  during  the 
previous  trajectory.  This  simple  restart  facility  thus  allows  the  contin- 
uation of  a trajectory  which  has  been  interrupted  by  a system  failure  or  for 
some  other  reason. 

If  the  trajectory  is  in  fact  new,  the  operator  is  asked  two  more  questions: 

(i)  MAIN  asks  whether  the  store  under  consideration  is  MK82  or  Karinga. 

This  conditions  the  system  to  take  data  from  the  appropriate  data 
store. 

(ii)  The  operator  is  then  asked  to  enter  a file  name,  which  specifies  the 
name  of  the  Constants  Block  set  up  by  ECCON,  as  described  in 
Section  4.1. 

MAIN  then  starts  the  "flight".  Subroutine  DAUX  is  called  to  initialize 
derivatives  after  which  the  dive  angle  is  calculated  by  means  of 
equation  (24)  and  the  aircraft  orientation  relative  to  earth  axes  is  determined 
using  equations  (28) . 

At  this  point  in  the  calculation  of  each  trajectory  point,  MAIN  tests  the 
parameter  IFREE,  which  is  set  positive  when  the  traverse  rig  console  button 
is  pressed  to  signify  a transition  to  a "bomb-alone"  trajectory.  If  IFREE 
has  become  positive,  a once-only  call  is  made  to  subroutine  TRFREE,  which 
converts  bomb  position  from  aircraft  axes  to  earth  axes,  the  latter  being 
more  meaningful  to  the  calculation  of  trajectory  to  impact. 

For  all  trajectory  points,  MAIN  then  calculates,  for  the  current  height, 
the  following  quantities: 

p = 10335.11  g(l  - 2.2559  x 10"5  h)  5-256103> 
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where  g 
h 

P 

T 


is  gravitational  acceleration  (m/s2.), 
is  height  (m) , 

is  static  pressure  (newton/m2), 
is  ambient  temperature  (°C) , 
is  velocity  of  sound  (m/s). 


is  free-stream  velocity  (m/s), 
is  free-stream  Mach  number, 
q is  dynamic  pressure  (newton/m2 ) and 
7 is  the  ratio  of  specific  heats. 

After  expressing  the  gravitational  force  on  the  bomb  in  aircraft  axes, 
MAIN  then  takes  one  of  twocoursesof  action,  depending  on  whether  the  bomb- 
alone  phase  is  in  operation. 


(a)  If  the  store  is  still  in  the  vicinity  of  the  aircraft, 

(i)  the  flow  properties  along  the  centre-line  of  the  aircraft, 
are  determined  (CALL  FLOW), 

(ii)  the  effective  incidence  and  sideslip  are  c.  lculated  using 
the  methods  of  Section  3.1  (CALL  CORECT) , 

(iii)  the  force  and  moment  coefficients  due  to  uniform  flow  are 
interpolated  from  the  disk  data  store  (CALL  CRLF), 

(iv)  the  increments  to  the  force  and  moment  coefficients  arising 
from  the  non-uniform  flow  are  calculated  (CALL  CINC),  and 
final ly, 

(v)  the  forces  and  moments  on  the  store  from  all  sources  are 
determined,  the  forces  being  in  aircraft  axes  and  the 
moments  being  in  bomb  axes. 

(b)  If  the  store  has  entered  the  "bomb-alone"  phase  of  its  trajectory, 

the  wind  velocity  components  in  bomb  axes  are  just  the  negatives  of  the 
bomb  velocity  components  in  earth  axes.  The  values  are  then  used 
in  steps  (a)  (iii)  and  (a)  (v)  above  to  determine  the  forces  and 
moments  on  the  store. 


Thus,  to  calculate  the  next  trajectory  point,  it  only  remains  to  invoke  a 
fourth-order  Runge  Kutta  integration  and  print  out  the  result. 

However,  the  steps  described  in  (a)  and  (b)  above  are  carried  out  only 
once  at  each  trajectory  point,  whereas  the  Runge  Kutta  process  requires  values 
of  the  forces  and  moments  at  several  places  within  the  time  interval  under 
consideration.  If  this  condition  is  not  met,  an  oversimplification  of  the 
integration  process  occurs  and  an  instability  is  produced  in  the  solution 
for  roll,  pitch  and  yaw. 

It  would  be  quite  impractical  to  calculate  forces  and  moments  at  every 
point  in  the  time  interval  required  by  the  Runge  Kutta  process,  because  of 
the  large  amounts  of  processing  time  involved.  Thus,  to  overcome  the 
instability,  a time-wise  parabola  is  fitted  to  the  three  most  recent  values 
for  each  of  the  three  forces  and  each  of  the  three  moments.  The  Runge  Kutta 
process  then  uses  each  parabola  (in  Subroutine  DAUX)  to  extrapolate  into  the 
current  time  interval  and  thus  to  produce  a more  accurate  estimate  of  each 
force  and  moment. 

The  process  described  above  is  repeated  to  calculate  each  trajectory 
point,  the  calculations  continuing  until  either  the  maximum  time  has  been 
exceeded  or  the  store  has  struck  the  ground. 

Specific  details  of  the  subroutines  used  by  the  Captive  Trajectory  Yawmeter 
System  are  given  in  subsequent  sections.  A listing  of  the  Main  Program 
appears  in  Appendix  II. 
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4.3  Subroutines  TRSET  and  TRFREE 

TRSET  is  a small  routine,  coded  in  Assembly  Language  (MACRO-15),  which 
initiates  the  Skip  Chain  (see  reference  11)  to  receive  interrupts  from  the 
traverse  rig  console  button.  It  is  called  by  MAIN  when  MAIN  is  first  loaded 
and  sets  the  parameter  IFREE  to  zero. 

When  the  traverse  rig  console  button  is  pressed,  control  is  transferred  to 
another  part  of  TRSET,  which  clears  the  flag,  turns  on  the  Program  Interrupt, 
sets  IFREE  to  +1  and  then  returns  control  to  the  interrupt  address. 

At  the  beginning  of  each  trajectory  point,  MAIN  tests  IFREE  to  see  if  it 
has  been  set  positive  (i.e.,  to  +1).  If  it  has,  then  MAIN  calls  Subroutine 
TRFREE  in  order  to  perform  the  transformation  to  a "bomb-alone"  trajectory. 

TRFREE  first  sets  the  parameter  FREE  to  +1.  It  then  performs  the 
transformation  of  coordinates  as  follows. 

If  a,  Ij  and  are  vectors  fixed  in  aircraft,  bomb  and  earth  axes  respect- 
ively, then 

a = CA  e (83) 


> - Si 


b = C a , (85) 

where  C^,  Cfi  and  C are  direction-cosine  matrices  as  defined  in  equation  (4). 
Equations  (83),  (84)  and  (85)  then  give 

CB  = C CA.  (86) 

Since  the  orientation  of  the  bomb  must  now  be  referred  to  earth  axes 
instead  of  aircraft  axes,  the  matrix  C is  replaced  by  the  matrix  Cg,  as 

defined  in  equation  (86) . 

Similarly,  if  x = (x,y,z)  is  the  vector  which,  until  the  transition,  has 
represented  the  position  of  the  bomb  centre  of  gravity  relative  to  earth 
axes,  then  this  can  be  expressed  as 

x = C^1  x + (87) 

where  x^  is  the  vector  which  represents  the  position  of  the  aircraft  centre 

of  gravity  relative  to  earth  axes.  The  components  of  x^  are  given  in 

equations  (25)  and  (26).  The  vector  x thence  describes  the  position  of  the 
bomb  centre  of  gravity  relative  to  earth  axes. 

Differentiation  of  equation  (87)  with  respect  to  time  gives 


C/1  u + 
A — 


V cos  (DIVE) 


-V  sin  (DIVE) 


Now  u signifies  the  velocity  of  the  bomb  centre  of  gravity  relative  to  earth 
axes  rather  than  relative  to  aircraft  axes. 
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Finally,  TRFREE  calls  Subroutine  DAUX  (to  set  up  the  new  derivatives), 
sets  the  gravitational  acceleration  components  (which  are  then  constant)  to 
earth  axes,  and  then  prints  a message  on  the  line  printer  to  signify  that  the 
"bomb-alone"  trajectory  has  begun. 

Listings  of  TRSET  and  TRFREE  are  given  in  Appendix  III. 

4.4  General  purpose  routines  (ARCTAN,  TRANS,  CTRANS  and  PARAB) 

Functions  ARCTAN  (Y,X)  calculates  the  angle  0 whose  sine  is  given  by  Y 
and  whose  cosine  is  given  by  X.  0 can  range  from  0 to  27T. 

TRANS  and  CTRANS  are  entry  points  to  a single  FORTRAN  routine.  For  given 
inputs  ROLL,  PITCH,  YAW  (Euler  angles),  INV  and  x,  TRANS  performs  one  of  the 
coordinate  transformations 


1*. 

II 

n 

IX 

(INV  = 

+ 1) 

(89(a)) 

x1  = C 1 x 

(INV  = 

-1), 

(89(b)) 

where  C is  calculated  according  to  equation  (4)  and  is  the  direction-cosine 
(orientation)  matrix  associated  with  the  transformation. 

CTRANS  is  a special  case  of  TRANS  in  that  it  transforms  specifically  from 
aircraft  to  bomb  axes  (INV  = +1)  or  vice  versa  (INV  = -1).  Since  the  nine 
elements  of  C in  this  case  are  all  dependent  variables  of  the  time  integ- 
ration, these  elements  are  always  available  for  use  in  equations  (89)  and 
thus  do  not  have  to  be  calculated. 

Subroutine  PARAB  simply  fits  a parabola 

f(x)  = a + b(x  - xA)  + c(x  - xA)2 

to  three  values  of  the  function  f(x)  and  stores  the  coefficients  in  the 
array  ABC. 

A listing  of  Subroutines  ARCTAN,  TRANS,  CTRANS  and  PARAB  is  given  in 
Appendix  IV. 

4.5  The  flow  measurement  and  probe  traverse  routines 

Subroutine  FLOW,  together  with  its  associated  subroutines  CALC,  TRAV, 
TRWAIT,  ADRD  and  PROBE,  is  responsible  for  finding  the  flow  properties  in 
bomb  axes  at  the  twelve  stations  along  the  store  centre-line. 

FLOW  first  calls  Subroutine  CALC  which 

(i)  transforms  the  coordinates  of  the  first  station  from  bomb  axes  to 
aircraft  axes, 

(ii)  Converts  the  coordinates  from  metres  to  inches, 

(iii)  divides  by  the  scale  factor  (SCALE)  to  convert  from  full  scale  to 
model  scale  length  measurements, 

(iv)  transforms  from  aircraft  axes  to  traverse  rig  axes,  and 

(v)  calls  Subroutine  TRAV  to  move  the  yawmeter  probe  to  the  position 
calculated. 

TRAV  is  an  Assembly  Language  routine  which  drives  the  traverse  rig  to  the 
wind  tunnel  coordinates  (X,Y,Z)  in  the  minimum  time.  It  uses  one  hardware 

register  to  drive  three  coordinate  axes  simultaneously  thus  making  its  opera- 
tion somewhat  complicated.  A full  description  of  TRAV  can  be  found  in 
Section  3.8.2  of  reference  5. 


When  the  probe  has  arrived  at  the  desired  position,  a probe  settling 
delay  of  six  seconds  is  allowed  (by  TRWAIT)  and  then  the  necessary  pressure 
measurements  are  taken  (by  ADRD).  Manifold,  pitot  and  the  two  differential 


WSRL-0005-TR 


24 


pressures  (DP13  and  DP24)  are  acquired  through  a Raytheon  Multiverter 
(multiplexed  analog-to-digital  converter),  while  the  total  pressure  is 
acquired  from  a digital  read-out  of  the  SI  wind  tunnel  Boulton-Paul  manometers. 

Having  acquired  the  pressure  data,  FLOW  again  calls  CALC  in  order  to 
start  the  probe  moving  towards  the  next  measurement  station.  Thus,  the 
probe  movement  to  the  next  station  and  the  processing  of  the  measurements  of 
the  current  station  can  proceed  in  parallel,  so  that  delays  can  be  kept  to 
a minimum. 

To  calculate  the  flow  properties  at  the  current  station  from  the  pressure 
measurements,  FLOW  calls  Subroutine  PROBE,  which  calculates  the  ratio  of 
pitot  pressure  to  manifold  pressure  and  prints  an  error  message  if  this 
ratio  is  less  than  one. 

PROBE  uses  the  pressure  measurements  from  the  yawmeter  probe  to  calculate 
a downwash  and  sidewash.  This  is  accomplished  by  means  of  three  probe 
calibration  data  stores  which  reside  on  disk.  The  three  files  are  as 
follows : 

(a)  TMACHT  PRB  resides  on  logical  unit  14  (.DAT  168)  and  is  a table  of 
Mach  number  as  a function  of  pitch  (PITCH)  and  the  pitot-to-manifold 
pressure  ratio  (PRATIO) , PITCH  ranging  from  0°  to  29°  in  steps  of  1° 
and  PRATIO  ranging  from  1 to  4.2  in  steps  of  0.05.  In  the  tables, 
PRATIO  ranges  more  rapidly  than  PITCH;  that  is,  the  first  65  table 
entries  are  for  a pitch  of  0°,  the  next  for  a pitch  of  1°,  and  so  on 
up  to  the  last  65  values  which  are  for  a pitch  of  29  . In  general, 
the  record  number  corresponding  to  a particular  PITCH  and  PRATIO  is 
given  by 


R = 65. PITCH  + 20(PRATI0-1)+1. 


(b)  PRATIO  PRB  is  on  logical  unit  3 and  is  a table  of  pitot-to-total 
pressure  ratio  (RATIO)  as  a function  of  pitch  ranging  from  0U  to  29 
in  steps  of  1°  and  Mach  number  ranging  from  0.4  to  1.4  in  steps  of 
0.1.  Again,  Mach  number  ranges  more  rapidly  than  pitch. 

In  general,  the  record  number  corresponding  to  a particular  PITCH  and 
MACH  number  is  given  by 


R = 11. PITCH  + 10. MACH-3. 


(c)  MACHCX  PRB  resides  on  logical  unit  12  (.DAT  148)  and  is  a table  of 
downwash  angle  and  sidewash  angle  as  functions  of  the  differential 
pressure  coefficients  (DC13,  DC24)  and  Mach  number.  DC13  and  DC24 
each  range  from  -2  to  +2  in  steps  of  0.1  and  Mach  number  ranges  from 
0.4  to  1.4  in  steps  of  0.1. 

Mach  number  ranges  most  slowly,  then  DC24  and  finally  DC13  ranges 
most  rapidly.  Thus,  the  first  41  entries  each  give  a downwash  and 
sidewash  for  a Mach  number  of  0.4  and  a DC24  of  -2,  the  next  41  entries 
correspond  to  a Mach  number  of  0.4  and  a DC24  of  -1.9,  and  so  on  up 
to  the  last  41  entries,  which  correspond  to  a Mach  number  of  1.4  and 
a DC 24  of  +2. 

In  general,  the  record  number  corresponding  to  a particular  DC13, 

DC24  and  MACH  number  is  given  by 


+ 


R 


412  (10. MACH-3) 


41l  10(DC24+2) ) + 10(DC23+2)+l . 
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Returning  to  the  operation  of  Subroutine  PROBE,  a pitot-to-manxfold 
pressure  ratio  is  calculated,  a pitch  equal  to  the  last  pitch  calculated 
(obviously  zero  for  the  first  time  through)  is  assumed,  and  the  file  TMACHT 
(as  described  in  (a)  above)  is  accessed  to  give  a corresponding  Mach  number 
(M)  . 

After  checking  that  the  Mach  number  is  in  the  range  0.4  to  1.4,  PROBE 
then  accesses  PRATIO  PRB  (as  in  (b)  above)  to  give  a pitot-to-total 
pressure  ratio  for  the  assumed  pitch  and  the  now  known  value  of  Mach  number. 
This  ratio,  together  with  the  measured  pitot  pressure,  then  gives  an 
estimate  of  total  pressure,  from  which  the  static  pressure  is  calculated 
from  the  formula 


P 


static 


Ptotal/(1  + C7  ' 1)M2/2)T/(7‘1) 


Dynamic  pressure  (q)  is  then  given  by  the  relation 


q 


T P 


static 


M2  . 


The  calculated  value  of  q,  together  with  the  measured  differential 
pressures  DP13  and  DP24,  enables  the  differential  pressure  coefficients  to  be 
calculated  from  the  formulae 


DC 13  = DP13/q 

and 

DC24  = DP24/q. 

The  file  MACHCX  PRB  (as  in  (c)  above)  is  then  used  to  give  values  of 
downwash  angle  and  sidewash  angle  as  functions  of  DC13,  DC24  and  M. 

The  downwash  and  sidewash  then  enable  a new  estimate  of  pitch  to  be 
calculated  and  the  process  is  then  repeated  using  this  updated  value  of  pitch 
The  iteration  is  said  to  converge  when  successive  approximations  for  pitch 
lie  within  0.15  of  each  other  but,  in  any  event,  an  error  message  will  be 
printed  and  execution  terminated  if  more  than  5 iterations  are  required. 

All  disk  accesses  are  performed  by  subroutine  READIN,  which  reads  from 
the  tables  entries  on  either  side  of  given  values  of  the  two  independent 
variables  and  then  forms  a two-dimensional  linear  interpolation. 

As  with  the  data  bank  of  store  coefficients,  the  accesses  to  the  disk 
to  obtain  probe  calibrations  are  carried  out  by  means  of  Direct-Access 
Input-Output  commands,  which  permit  the  user  to  directly  reference  any  record 
in  a file  without  indexing  from  the  file's  beginning  up  to  the  desired  record 
Therefore,  all  calls  to  READIN  must  include,  not  only  the  logical  unit 
number  to  define  the  file  to  be  accessed,  but  also  a record  number  which  is 
calculated  from  the  values  of  the  independent  variables  under  consideration. 
PROBE  thus  contains  the  necessary  algorithms  for  computing  appropriate 
record  numbers  associated  with  the  files  TMACHT,  PRATIO  and  MACHCX. 

When  PROBE  has  calculated,  relative  to  probe  axes,  the  downwash  angle 
(DOWNP)  and  sidewash  angle  (SIDEP)  for  the  current  store  measurement  station, 
control  is  then  returned  to  Subroutine  FLOW. 

FLOW  then  uses  the  Mach  number  calculated  by  PROBE  to  determine  the  local 
velocity  and  then  the  components  of  local  velocity  relative  to  probe  axes. 
These  are  given  by 
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N aoo)'  = (1  + — 

-jj—  Hx>2  )/(l  + M2) 

(90) 

1 VI 

= Ma 

(91) 

Up  = - 1 VI  / yj  1 + 

tan2 (DOWNP)  + tan2  (SIDEP) 

(92(a)) 

vp  = Up  tan(SIDEP) 

(92(b)) 

Wp  = Up  tan  (DOWNP),  (92(c)) 


where  a is  the  local  velocity  of  sound  and  the  subscript  °°  refers  to  free- 
stream  values. 

V is  the  wind  vector  and  is  shown  in  figure  5 together  with  the  probe  axes 
and  the  corresponding  velocity  components  Up,  Vp  and  Wp. 

The  vector  (Up,  Vp,  Wp)  is  transformed,  through  the  Euler  angles  probe 

roll  and  probe  pitch,  to  wind  tunnel  axes,  then  to  aircraft  axes  and  finally 
to  bomb  axes.  The  final  velocity  components  (UMS,  VMS,  WMS)  returned  by 
FLOW  therefore  represent  the  components  of  wind  velocity  relative  to  bomb 
axes. 

FLOW  repeats  the  entire  process  for  each  of  the  twelve  measurement  stations. 
To  save  time,  the  measurement  stations  are  visited  in  the  sequence 


(1,2,3,  12)  for  odd  trajectory  points  and  in  the  sequence 

(12,11,  1)  for  even  trajectory  points  - a "zig-zag"  effect. 


Listings  of  Subroutines  FLOW,  CALC,  TRAV,  TRWAIT,  ADRD,  PROBE  and  READIN 
are  given  in  Appendix  V. 

4.6  The  uniform  flow  routines  CORECT  and  CREF 

Section  3.1  describes  how  uniform  flow  conditions  are  calculated  at  a 
reference  point. 

Subroutine  CORECT  uses  the  methods  of  Section  3.1  to  calculate  the 

effective  incidence  (a  ) and  the  effective  sideslip  (j3  ),  given  by 

0 0 


= a + 2y  (93) 

and 

PQ  = P * 27,  (94) 

where  u,  0 and  7 are,  respectively,  the  incidence,  sideslip  and  camber 
defined  in  Section  3.1. 

CORECT  then  returns  control  to  the  Main  Program,  which  calculates  total 
pitch  (0,p)  and  total  roll  (>^>T)  according  to  the  formulae 

0T  = tan'1  {yf  (V/U)2  + (W/U)2)  (95) 


and 


<0 


T 


ARCTAN  (V/U,  W/U). 


(96) 


r 
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Equations  (95)  and  (96)  derive  from  the  equations  (see  figure  4) 

U = -IVl  cos  eT  (97(a)) 

V = -I  V|  sin  0T  sin  ^ (97(b)) 

and 

W = -IVl  sin  0T  cos  (97(c)) 

where  V is  the  wind  vector  and  U,V,W  are  the  components  of  the  wind  vector 
relative  to  bomb  axes. 

The  next  step  is  to  call  subroutine  CREF,  which  uses  total  pitch  and 
total  roll  to  read  the  appropriate  force  and  moment  coefficients  from 
logical  unit  13  (.DAT  158)  of  the  disk. 

This  uniform  flow  data  store  is  named  'MK82CX  STC’  for  the  Mk82  Bomb 
data  and  is  named  1 KARCX  STC'  for  the  Karinga  store  data.  Each  file 
consists  of  8 x 13  x 16  = 1664  records,  each  record  containing  the  10 
coefficients  CX,  CY,  CZ,  CL,  CM,  CN,  CLP,  CMQ,  CNP  and  CYP  for  given  Mach 
number.  Total  Pitch  and  Total  Roll.  The  last  four  coefficients  are 
dynamic  derivatives  and  are  given  by 


CLP  = 

2V  6CL 

D 3p 

(Roll  Damping) 

(98(a)) 

- 

CMQ  = 

2V  3CM 

D 9q 

(Pitch  Damping) 

(98(b)) 

CNP  = 

2 v acN 

D 3p 

(Magnus  Moment) 

(98(c)) 

and 

CYP  = 

2V  3CY 

D 3p 

(Magnus  Force) 

(98(d)) 

where  D is  maximum  cross-sectional  diameter,  V is  local  velocity  and  p,  q,  r 
are  bomb  Euler  rates  relative  to  total  axes. 

The  range  of  Mach  numbers  is  0.4,  0.5,  0.6,  0.7,  0.8,  0.85,  0.9  and  0.95 
(this  range  is  currently  being  extended  to  Mach  1.2).  Total  roll  ranges  from 
-45  to  +45  in  steps  of  7.5  while  total  pitch  ranges  from  0°  to  30°  in  steps 
of  2 . 

Mach  number  ranges  most  slowly,  then  total  roll,  and  finally  total  pitch 
ranges  the  most  rapidly.  Thus,  the  first  16  entries  correspond  to  a Mach 
number  of  0.4  and  a total  roll  of  -45°,  the  next  16  to  M = 0.4  and  total 
roll  = -37.5  , and  so  on  up  to  the  last  16  values  which  correspond  to  a Mach 
number  of  0.9  and  a total  roll  of  +45°. 

In  general,  the  record  number  corresponding  to  a particular  Mach  number, 
total  roll  (<p  ) and  total  pitch  (0^,)  is  given  by 


f, 


! 


! 


= 208  (M  - 1)  + 16[  0PT  + 45)  / 7 . 5]  + [0.J./2]  + 1, 


R 
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where  l ] denotes  the  integer  part  and  M is  defined  by 

M = | 10. MACH  - 3)  for  0.4  < MACH  <0.8 


and 


M = [ 20. MACH  - 1 1 J for  0.8  < MACH  < 0.95. 


CREF  calls  READIN  twice  to  determine  the  flow  values  at  the  Mach  number 
table  entries  above  and  below  the  required  Mach  number.  It  then  interpolates 
over  Mach  number  to  give  the  final  values  of  the  coefficients.  The  total 
interpolation  therefore  consists  of  two  two-dimensional  interpolations  over 
pitch  and  roll  (one  for  each  call  to  READIN)  plus  a one-dimensional  inter- 
polation over  Mach  number. 

Allowing  for  the  dynamic  derivatives,  CREF  then  calculates  the  coefficients 
relative  to  total  axes  as 


and 


TCX  = CX 

(99(a)) 

TCY  = CY  + 

3CY 

^ dp 

(99(b)) 

TCZ  = CZ 

(99(c)) 

TCL  = CL  + 

3CL 

P 3p 

(99(d)) 

TCM  = CM  + 

3 CM 
q 3q 

(99(e)) 

TCN 


CN  + 


dCN 
P dp 


3CN 
r 3r  ’ 


(99(f)) 


where  CX,  CY  etc.  are  read  from  the  data  store  and  the  dynamic  derivatives 
are  obtained  via  equations  (98)  also  from  the  data  store.  It  is  assumed 

that  = |^.  Finally,  CREF  transforms  the  six  calculated  coefficients 

from  total  axes  to  bomb  axes  and  then  returns  control  to  the  Main  Program. 
Listings  of  CORECT  and  CREF  are  given  in  Appendix  VII. 

4.7  The  non-uniform  flow  routine  C INC 


Subroutine  CINC  uses  the  methods  of  Section  3.2  to  calculate  non-uniform 
flow  increments  to  the  coefficients  calculated  in  Section  4.6. 

Since  Section  3.2  fully  discusses  the  techniques  involved,  it  is  sufficient 
to  refer  the  reader  to  that  section  rather  than  repeat  the  description  here. 

The  derivatives  ^ required  by  equations  (76)  and  (81)  are  calcul- 

ated in  the  following  way: 
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(i)  If  the  derivative  is  required  at  the  point  x^,  then  CINC  fits 
the  parabola 


F = a + bx  + cx2 


successively  through  the  points  (x^  2 , xA  ^ , x^),  x^  x^+1) 

and  (x.,  x.+1,  x.+2). 

(ii)  Each  estimate  of  the  derivative  is  then  given  by 


dF 

dx 


b + 2cx. 


(iii)  The  arithmetic  mean  of  the  three  estimates  is  then  calculated  to 

dF 

give  the  final  estimate  of  ^ at  the  point  x^.  This  averaging 

technique  gives  good  results  and  seems  to  obviate  the  problems  which 
arise  in  a simple  curve  fit  when  the  points  are  unevenly  spaced. 

When  CINC  has  calculated  the  non-uniform  flow  force  and  moment  increments 
at  each  of  the  bomb  measurement  stations  (x^),  it  integrates  each  contribution 

with  respect  to  x to  give  the  non-uniform  flow  and  moment  increments  for  the 
entire  bomb. 

A listing  of  Subroutine  CINC  appears  in  Appendix  VII. 

4.8  The  integration  routines  INTM  and  DAUX 

Subroutine  INTM  is  simply  a classical  fourth-order  Runge-Kutta  integrator 
similar  to  those  which  have  been  in  existence  for  some  years  now. 

Briefly,  let  the  system  of  equations  to  be  solved  be  given  in  the  form 


x = f(x,t). 


If  the  point  t has  been  reached  and  x(t)  = x,  INTM  calculates  in  succession 
the  quantities 


ko 

= hf (x,  t). 

k, 

= hf  (x  + ’sko  , 

t + 4h), 

k2 

= hf(x  + 4k,  , 

t + 4h) , 

and 

k3  = hf (x  + k2 , t + h) , 


where  h is  the  time  step. 
Then 


x (t  + h) 


x(t)  + (Jq,  + 2k,  + 2k2  + k3)/6. 
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In  fact,  INTM  integrates  18  first-order  differential  equations  at  each 
time  step.  To  do  this  it  uses  the  first  75  elements  of  the  main  program 
matrix  T. 

T(l)  is  used  to  store  time,  T(2)  is  the  current  time  and  T(3)  is  the 
time  step. 

T(4)  to  T(21)  are  the  dependent  variables  to  be  integrated,  where 

(a)  T (4)  to  T (12)  are  the  elements  C (1, 1) , C(2,l),  C(3,l),  C(l,2),  

C(3,3)  of  the  direction-cosine  matrix  C as  defined  in  equations  (4) 
and  (47) . 

(b)  T(13)  to  T(15)  are  X,  Y,  Z,  the  position  of  the  bomb  centre  of 
gravity  relative  to  aircraft  axes, 

(c)  T(16)  to  T(18)  are  U,  V,  W,  the  velocity  of  the  bomb  centre  of  gravity 
relative  to  aircraft  axes, 

(d)  T(19)  to  T(21)  are  the  angular  momentum  components  of  the  bomb 
relative  to  earth  axes,  but  expressed  in  bomb  axes, 

(e)  T(22)  to  T(30)  are  the  rates  of  change  of  the  direction-cosines  in 

(a)  above,  as  given  by  equation  (54), 

(f)  T(31)  to  T(33)  are  the  rates  of  change  of  X,  Y,  Z in  (b)  above, 

(g)  T(34)  to  T(36)  are  the  rates  of  change  of  U,  V,  W in  (c)  above,  and 

(h)  T(37)  to  T(39)  are  the  rates  of  change  of  the  angular  momentum 

components  in  (d)  above. 

In  addition,  INTM  uses  T(40)  to  T(75)  as  working  storage. 

Whenever  INTM  wishes  to  calculate  the  derivatives  corresponding  to  the 
18  first-order  differential  equations  it  calls  Subroutine  DAUX. 

DAUX  first  updates  the  angular  velocities  p,  q,  r from  the  angular 
momentum  components  by  means  of  equation  (40). 

It  also  re-normalizes  the  elements  of  the  direction-cosine  matrix  by 
means  of  the  formula 


C1.  . 
ij 


(i  = 1,2,3), 


The  truth  of  this  relationship  follows  from  equation  (4),  which  shows  that 
the  sums  of  the  squares  of  the  elements  in  any  row  or  column  of  the 
direction-cosine  matrix  C is  unity. 

Next  DAUX  updates  the  bomb  and  aircraft  rotation  matrices  and  A 

from  the  angular  rates  p,  q,  r and  p q , r respectively,  according  to 

t\  J\ 

equation  (15). 

Then  the  following  derivatives  are  calculated: 

(i)  The  direction-cosine  matrix  derivatives  (C) , one  for  each  of  the 
nine  elements,  given  by  equation  (54)  as 


= fl.c-csi.. 

B A 


(100) 


(ii)  The  derivative  ()T)  of  (X,Y,Z),  given  by  equation  (34)  as 
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(iii)  The  derivative  (u)  of 


u = F^/m  - 


(iv)  The  derivatives  Qu)  of  the  angular  momentum  components,  given  by 
equation  (46)  as 

te  - n b V (103) 

As  noted  in  Section  4.2,  the  forces  (F^)  in  equation  (102)  and  the 

moments  (Mg)  in  equation  (103)  are  calculated  by  DAUX  via  a parabolic 

extrapolation  into  the  current  time  interval,  in  order  to  overcome  an 
inherent  instability  in  the  integration  process. 

It  should  be  noticed  that,  once  the  "bomb-alone  trajectory"  phase  has 
begun,  the  variables  describing  the  state  of  the  bomb  are  referred  to  earth 
axes,  and  therefore  DAUX  then  ignores  the  final  terms  of  equations  (100), 
(101)  and  (102).  That  is,  the  3x3  matrix  £2  A = 0. 

A listing  of  Subroutines  INTM  and  DAUX  is  given  in  Appendix  VIII. 


(U,  V,  W) , given  by  equation  (39)  as 


Mg  sin (ATTACK)' 
0 

■Mg  cos  (ATTACK)/ 


+ , u . 

A — 


(102) 


5.  CONCLUDING  COMMENTS 

The  Captive  Trajectory  Yawmeter  System  appears  to  be  a successful  means  of 
simulating  store  separation  from  aircraft.  Results  to  date  are  most  encouraging. 
The  technique  of  using  a yawmeter  probe  in  place  of  a model  mounted  upon  a sting 
decreases  problems,  such  as  scale  effects,  associated  with  small  wind  tunnels. 

Since  a typical  trajectory  requires  approximately  30  min  of  running  time,  the 
CTYS  is  aimed  at  small  continuous  flow  facilities.  Because  a settling  delay  of 
approximately  5 s is  required  at  each  of  the  12  measurement  stations,  probe 
settling  time  therefore  accounts  for  one  minute  per  trajectory  point.  When,  in 
addition,  the  traversing  time  of  the  rig  is  taken  into  account,  it  is  obvious 
that  the  quoted  total  time  of  30  min  cannot  be  reduced  for  any  reasonable  traject- 
ory, despite  the  fact  that  data  store  accesses  have  been  streamlined  to  the  utmost 
and  are  carried  out  in  parallel  with  the  traversing  and  settling  operations. 

Thus,  the  CTYS  would  be  clearly  unsuited  to  blowdown  tunnels  and  very  expensive 
to  use  in  large  continuous  flow  facilities. 

Two  extensions  of  the  system  are  currently  underway.  Firstly,  a model  of 
the  ejectors  is  being  developed (ref .4)  so  that  store  position,  attitude  and 
rates  at  the  end  of  the  ejector  stroke  can  be  fed  automatically  to  the  program. 

At  present,  these  must  be  supplied  by  the  experimenter  as  input  parameters. 

Secondly,  it  has  been  found(ref.2)  that  the  yawmeter  probe  can  give  spurious 
readings  when  it  is  positioned  too  close  to  any  part  of  the  aircraft  surface, 
such  as  bomb  racks.  Consequently,  off-line  measurements  are  currently  being 
made,  to  provide  a data  store  of  aerodynamic  coefficients  as  a function  of 
position  and  attitude  whilst  the  probe  is  located  within  this  difficult  region. 
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NOTATION 

AROLL, APITCH.AYAW 

Euler  Roll,  Pitch  and  Yaw  of  aircraft  relative  to  Earth 
axes 

C 

direction-cosine  (coordinate  transformation)  matrix 

°C 

cross-flow  drag  coefficient 

C.  . 

1J 

ij^*1  element  of  C 

D 

cross-sectional  diameter  of  bomb  (metres) 

DIVE 

aircraft  dive  angle 

F 

external  forces  on  bomb 

IXX,IYY’IZZ 

moments  of  inertia  of  bomb  about  its  principal  axes 
(kgm  m2) 

K 

mean  chord  of  aerofoil 

M 

Mach  number 

M 

external  moments  on  bomb  (newton-m) 

S 

cross-sectional  area  of  bomb  (=  n D2/4) 

T 

ambient  temperature 

u.v.w 

velocity  components  of  bomb  C.G.  relative  to  Aircraft 
axes  (m/s) 

X,  Y , Z 

coordinate  axes;  position  of  bomb  C.G.  (metres) 

X*,Y*,Z* 

position  of  aircraft  C.G.  (metres) 

XF,YF,ZF 

coordinates  of  bomb's  reference  fin  relative  to 

Aircraft  axes  (metres) 

- 

XN.YN.ZN 

coordinates  of  bomb's  nose  relative  to  Aircraft  axes 

i 

Y,Z 

cross-force  components  (see  Section  3.2) 

a 

speed  of  sound  (m/s) 

h 

[ 

a 

vector  fixed  in  Aircraft  axes 

b 

vector  fixed  in  Bomb  axes 

e 

vector  fixed  in  Earth  axes 

f 

potential  cross-force  (see  equation  (75)) 

g 

gravitational  acceleration  (9.80665  m/s) 

h 

height  of  aircraft;  mean  height  of  aerofoil 

33 


WSRL-0005-TR 


h 

m 

P 


p»q.r 


q 

t 

u 


x,y,z 

£2 


a 


a 

e 


0 

0 


e 


7 


v*  ,o*  ,P* 


i 

Mg 

p 

00 

Subscripts 

A 

B 

E 

P 

T 

V 

u 

0 

oo 


angular  momentum  vector  (newton-m-s) 
mass  of  the  bomb  (kgm) 
static  pressure  (newton/m2) 

angular  rates  of  rotating  system  about  its  own  X,Y,Z 
axes  respectively 

dynamic  pressure  (h  P V2 ) newton/m2 ) 
time  (s) 

see  equation  (33) 

position  of  aircraft  C.G.  relative  to  earth  axes 

angular  rotation  matrix 

incidence 

effective  incidence  (u  + 27) 
sideslip 

effective  sideslip  Q3  + 27) 

wing  camber  (=h/c);  ratio  of  specific  heats 

Euler  Roll,  Pitch  and  Yaw 

Projected  Roll,  Pitch  and  Yaw 

crossflow-drag  proportionality  factor  (see  Section  3.2) 
centripetal  acceleration  of  aircraft  trajectory 
local  air  density;  radius  of  curvature 
angular  rotation  rate  vector 

Aircraft  axes 
Bomb  axes ; Buoyancy 
Earth  axes 
Potential;  Probe 
Total  axes 
Viscous 
Uniform  Flow 
initial  conditions 
free-stream  conditions 
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APPENDIX  I 

LISTINGS  OF  THE  CONSTANTS  EDITOR  (EDCON) 


Li 


C EDCON  DATE;  2*-6-77 

C 

C CTS  CONSTANTS  EDITOR 

CALLING  SEQUENCE: 

t A RK  A -5 
♦GLOAD 

LOADER  V3ABBB 
> -EDC  ON  <ALTMODE  > 

r 

DIMENSION  F ILE<  2 > 

DIMENSION  XMS<  12  >,  OM$( 12  >,  AMS<  12) 

DIMENSION  T(32B>, IT<5  >, TS(  17 > 

DIMENSION  C<3 , 3 > , P V( 3 )/ DELPV<  3 > , PAV( 3) 

DIMENSION  FSV(3,6),TIMES(3>,ABC(3,6> 

C 

C*****TINE,  INTEGRATION  TIME  STEP 

EQUIVALENCE  (T(2),TIME),(T(3>,DT> 

C***»*D  IRECT  ION  COSINES  DEFINING  ORIENTATION  OF  BOMB  WITH  RESPECT  TO 
C AIRCRAFT. 

EQUIVALENCE  ( T < 4 ) , C 1 1 , C ) . ( T( 5 > , C2 1 > , ( T ( S ) , C3 1 > 

EQUIVALENCE  ( T ( 7 ) , C 1 2 ) / < T<  8 ) , C22  > , < T ( 9 ) , C 32  ) 

EQUIVALENCE  (T(1B),C13>,(T(11>,C2  3>,(T(12>,C33) 

C *****  p OS  I T I ON  OF  BOMB  CENTRE  OF  GRAVITY  IN  AIRCRAFT  AXES 
EQUIVALENCE  (T(13>,X),(T(14),Y),(T(13>,Z> 

C ***** V EL OC I TY  OF  BOMB  C.G  IN  AIRCRAFT  AXES. 

EQUIVALENCE  (T(16>,U>,(T(17>,V>,(T(18>,W> 

C * **  ** ANGUL  A R MOMENTUM  OF  BOMB  IN  BOMB  AXES 

EQUI  VALENCE  < T(  I 9 > , NX  ),  ( T<  28  > , H Y ) , < T(  2 1 ),  HZ  ) 

C*****RATE  OF  CHANGE  OF  DIRECTION  COSINES  ABOVE 

EQUIVALENCE  ( T ( 22  ) , DC  1 1 ) , ( T<  2 3 ) , D C 2 1 ) , < T<  24  > , D C3 1 ) 
EQUIVALENCE  ( T< 2 5 ) , DC  1 2 ) , < T<  26  > , D C 22  ) , < T<  27  ) , D C3 2 ) 
EQUIVALENCE  < T ( 2 8 ) , DC  1 3 ) , < TC  2 9 ) . D C 23  > , ( T(  3B  ) , D C3 3 > 

C*****RATE  OF  CHANGE  OF  X,Y,Z, ABOVE 

EQUIVALENCE  ( T ( 3 1 ) , DX  ) , < TC  32  ) , D Y ) , < T < 3 3 ) , DZ  ) 

C *****  R ATE  OF  CHANGE  OF  U,V,W  ABOVE 

EQUIVALENCE  ( T ( 3 4 > , DU  ) , ( T<  35  ) , D V ) , C T ( 3 6 ) , DW  > 

C*****RATE  OF  CHANGE  OF  HX,HY,HZ  ABOVE 

EQUIVALENCE  < T < 3 7 ) , DH X ) , < T < 3 8 ) . DH Y ) , < T < 39  ) , D HZ  > 

C*****T<40)  TO  T( 75  > RESERVED  FOR  IN  T M 

,;*****A  IRCRAFT  ROTATION  ANGLES  W.R.T  EARTH  AXES 

EQUI VALENCE  ( T<  7 6 > , AROL L ) , ( T < 77  ) , A P I TCH  ) , < T ( 78  ) , A YA U > 
C*****BOMB  MOMENTS  OF  INERTIA 

EQUIVALENCE  <T(?9>,AXX),(T(8Bi.AYY>,(TC81),AZZ> 

C*****BOMB  MASS  (KGM) 

EQUIVALENCE  ( T ( 8 2 > , ST  MA  S S > 

C*****PROBE  ROLL  AND  PITCH  ANCLES  W.R.T.  TRAVERSE  AXES 
EQUIVALENCE  ( T< 83  > , APTROL  > , < TC 84  ) , APTP I T ) 

C *****  T IME  LIMIT  FOR  DROP, AIR  TEMPERATURE  (DEG  C>, 
t*****BOMB  C.G.  <. METRES  FROM  T A I L > , A I RC RA F T HE  I GH  T(  M E T RE  S > , 

C»**** BOMB  SCALE  (AS  IN  I / S CA LE  ) . T ( 9 B > SPARE, SPEED  OF  SOUND, 
C*****ACCELERATI ON  DUE  TO  GRAVITY, T( 93)  SPARE 

EQUIVALENCE  (T(85),TLIHIT>,(T(8G),TEMP),(T(87>,XCG>, 

1 ( T( 88) , HGHTB >, ( T( 89 >. SCALE >, 

2 ( T C 9 1 >, VSOUND  > , ( T(92), GRAVAC) 

C*****TOTAL  PITCH  AND  ROLL  OF  80MB  IN  FLOWFIELD  AT  REFERENCE  POINT 
EQUIVALENCE  ( T ( 9 4 ) , THET 0 T ) , ( T ( 95  ) . PH  I T OT  > 


4 


c*****free  stream  dynamic  pressure 

EQU I VALENCE  ( T ( l 04  ) , Q DP  8 ) 

C*****WIND  FORCES  ON  BOMB  IN  AIRCRAFT  AXES. 

EQUIVALENCE  (T(105),MX),(T(1B6),MY),(T(1B7),MZ) 

C*****WIND  MOMENTS  ON  BOMB  IN  BOMB  AXES. 

EQUIVALENCE  <T(1B8).ML>.(T(109),MH),(T(11B),MN) 

C *****  ft  IRCRAFT  ROLL, PITCH, YAM  RATE  ABOUT  AIRCRAFT  AXES 

EQUIVALENCE  ( T ( 1 1 1 ) , P A , P A V ) , < T<  1 1 2 ) , QA  ) , < T<  1 1 3 ) , R A > 

C 

C***»*  FREE  = -1(4-1.  > IF  BOMB-ALONE  TRAJECTORY 
C IS  NOT  (IS)  OPERATING . 

C 

EQUIVALENCE  ( T< 1 14  ) . FREE  ) 

C*****ALPHA  AND  BETA  OF  BOMB  USED  FOR  CALCULATING  COEFFICIENTS 
EQUIVALENCE  <T(115),ALPHA>,(T(ii6),BETA) 

C***-**A  IRCRAFT  ACCELERATIONS  IN  AIRCRAFT  AXES. 

EQUIVALENCE  (T(U7),AX),(T(118).AY),(T(119>,AZ) 

C*****FORCES  ON  BOMB  C.G.  IN  AIRCRAFT  AXES. 

EQUIVALENCE  (T(12B),FX.,(T(121),FY>,<T(122>,FZ) 

C*****H0MENTS  ON  BOMB  ABOUT  B0M8  X,Y,Z  AXES 

EQUIVALENCE  <T(I23),FL>,(T(124),FM),(T(125),FN> 
C*****ROTATION  ANGLES  DEFINING  ORIENTATION  OF  B0M8  M.R.T  AIRCRAFT 
EQUIVALENCE  ( T ( 1 26  > , R OL L ) - ( T ( 1 2 7 ) , P I TC H ) , ( T < 1 2 8 ) . Y A W ) 
C*****R0LL, PITCH, YAW  RATES  OF  BOMB  ABOUT  BOMB  X,Y,Z  AXES 
EQUIVALENCE  (T(129>,P,PV>,(T(13B),Q>,(T(131),R> 

C*****PRES3URE  AT  H E I G H T , D Y N A M I C PRESSURE  AT  HEIGHT 
EQUIVALENCE  (T(132),PRESS),(T(133),QDPALT) 

C*****DEGREES  TO  RA D I A NS , R A D I AN S TO  DEGREES 

EQUIVALENCE  (T(134),DTR),(T(135),RTD) 
f**»** VELOCITY  OF  AIR  AT  -INFINITY- 
EQUIVALENCE  < T( 136 >.VELINF  ) 

C *****  A IRCRAFT  SELF -ROTAT I ON-RATE  MATRIX. 

EQUIVALENCE  <T(137),AR11),(T<138),AR21),<T(139),AR31) 
EQUIVALENCE  < T < 1 40  ) , A R 1 2 ) , < T ( 1 4 1 > ■ AR 22  ) , ( T(  1 42  ) , A R 3 2 ) 
EQUIVALENCE  <T(143),AR13),(T(144>,AR23),(T(145),AR3  3) 
C*****60MB  SELF-ROTATION-RATE  MATRIX. 

EQUIVALENCE  (T(146),BRli>,(T<147),BR21),(T(148),BR31> 
EQUIVALENCE  ( T < 1 49  ) , B R 1 2 > , < T < 1 SB  ) , BR 22  ) , < T<  1 5 1 ) , B R3 2 > 
EQUIVALENCE  ( T ( 1 52  ) , B R 1 3 > , ( T < 1 5 3 ) , BR 23  > , < T<  1 54  > , B R3 3 ) 

C*****T(  155  ) TO  T ( 1 6 0 > SPARE. 

C 

C*****GRAVI TY  FORCE  IN  AIRCRAFT  AXES. 

EQUIVALENCE  (T(16l),GX),<T(l62),GY),(T<163),GZ) 

C*****DIVE  ANCLE  OF  AIRCRAFT 

EQUIVALENCE  ( T < 1 64  ) , D I V E ) 

C*****pl,TWOPI, CONVERSION  FROM  METRES  TO  INCHES 

EQUIVALENCE  (T(165),PI>,(T(166),TM0PI),(T<167),CMI) 

C*»**  CENTRIPETAL  ACCELERATION  OF  AIRCRAFT  (IN  C'S) 

EQU I VALENCE  ( T(  1 68  ) , GF  ) 

C*****AT  ALL  MEASUREMENT  STATIONS  A RE  AS , D I AM E TE RS , D Y NA M I C PRESSURE, 
C*****NACH  NUMBER, U,V,W  (MIND  VELOCITY  CO M PO HE N TS , B 0 MB  AXES), 

C *****  AND  MEASUREMENT  STATIONS  (C  C =0.NOSE  + VE  , I E , BOMB  AXES) 
EQUIVALENCE  ( T ( 1 69  ) , AMS  ) , ( T(  l 8 1 ) , D MS  ) , ( T(  1 9 3 ) , QM S ) . 
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1 ( H2  05  ) - RMS  ) , ( T(  2 17  >,  UMS  )»  ( T<  22  9 >,  VMS  >,  < i ( 24  1 > , WMS), 

2 ( T ( 253  > , XMS  ) 

C *****  F RE  E STREAM  AERODYNAMIC  COEFFS  BOMB  AXES 

EQUIVALENCE  ( T ( 2 65  ) , CX ) , ( T ( 26 6 ) , C Y ) , < T < 26  7 ) , CZ  > , 

1 (T(2  68  >,CL>,(T(269>,CM>,(T(27B>,CN> 

C***«* INCREMENTS  TO  COEFFS  DUE  TO  NON-UNIFORM  FLOW 

EQUIVALENCE  < T < 2 7 1 ) , D CX  ) . ( T<  2 72  ) , D CY  > . < T<  27 3 > . DC Z ) . 

1 ( T(274  > ,DCL  >,(  T(  275  >,  DCM  >,  ( T(  276  >,  DCN  > 

c*****probe  pressure  measurements 

EQUIVALENCE  (H277),0P13)»(T(2?8>;DP24), 

1 ■ H 279  ) .PM  AN  ),  < 7<  280)  - PPIT  >,(  T<  2 8 i ),  T PRESS  ) 

C AREA  TO  SAVE  COEFFICIENTS  AT  PREVIOUS  2 TRAJ  POINTS 

EQUIVALENCE  ( T<  282  ) , FSV  ) , ( T(  300  ),  T IMES  > , < T<  3B3  ),  ABC  > 


* * * **  N 0 OF  MEASUREMENT  STATIONS, REF  PO INT , I NTEGRAT I ON  LIMITS. 
EQUIVALENCE  < I T<  1 ) , NM S ) , ( I T<  2 ) , NR E F ) , < I T<  3 ) , NB AC K > 
EQUIVALENCE  ( I T( 4 ) , NF RO N T ) , < I T< 5 ) , NPR I N T ) 


EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 

EQUIVALENCE 


(TSC  4), ROLLS), (TS(5), PIT  CHS).(TS(6>, YAWS) 

< TS<  7 >,  PS  ),  < TS(8  ),  QS),  < TS<  9 ),  RS  > 

( TS<  IB  ) , ATT  ACS  >,  < TS<  1 I ),  DI  VES  ) 

( TS(  12  >,PAS  ),(  TS(  13 ),  QAS ), < TS< 1 4 >, RAS  ) 

(TS(  15),SPTR0L),(TS<16),SPTPIT) 

( DELPV( 1 ),DELP ), < DELPVC2  >. DEL 8) , (DELPV(  3), DELR  ) 


CONSTANT  DATA: 

DATA  A NS  N0/5HN  / , ANS YES/5HY  /,NM$/12/ 

DATA  NREF/2/, N8ACK/3/ , NFRONT/ 12/ 

DATA  PI/3 . 1 4159265/ ,CMI/39 . 3 7B079/,CRA VAC/9  80665/ 
DATA  DTR/B. B1 7453293/, RTD/57 . 29578/ 

C 

TAN( X)=SIN<  X > / CO  S( X ) 


1 F0RMAT(///1X' CTS  CONSTANTS  E D I T OR ' // 1 X ' EN TE R FILE  NAME  ') 

2 FORMAH  A5  , A4  ) 

3 F0RMAK//1X, ' CREATING  NEW  FILE?') 

4 FORMAT(//lX'***BOMB  PARAMETERS*** ' ) 

714  F0RMAT(/1X' MOMENTS  OF  INERTIA') 

5 FORM AT< 1 X ' I X (KGM  M2):'Fi2.5) 

b F0RNAHF12. 6) 

7 FORMAT ( IX ’ I Y (KGM  M2):'F12.5) 

8 FOR  MAH  IX'  IZ  (KGM  M2)s'F12.5> 

9 FORMAK  / 1 X'  MASS  ( K G M > : ' F 1 2 . 5 > 

IB  FORM  AT ( / 1 X' C G POSN  (METRES  FROM  T A I L ) : ' F 1 2 . 5 ) 

11  F0RMAH/1X'  BOMB  SCALE  (S  AS  IN  1/S):'F12.5> 

12  F0RHAT(/1X'ENTER'. 14, ' MEASUREMENT  STATIONS  (M>'//> 

13  F0RMAT(//1X'ENTER' , 14,'  DIAMETERS  AT  HEAS.  STATIONS  (M)'/) 

14  FORMAT( // IX ' POSH  OF  BOMB  C.G.  (AIRCRAFT  AXES)') 

15  FORM AT( 1 X ' X ( M E T RE  S > s ' F 1 2 . 5 ) 

16  FORMATdX'Y  ( NET  RES  ) : ' F 1 2 . 5 > 

17  FORMAT ( IX'Z  ( N E T RE  S > : ' F 1 2 5 ) 

18  F0RMAT(//1X'  VELOCITY  OF  BOMB  C.G  (AIRCRAFT  AXES)') 

19  FOR  MAH  lX’U  (H/S):'FI2  5) 

2 B FORMAT ( 1 X ' V < M/S > : ' F 1 2 . 5 ) 

21  FORMAT(  lX'W  <M/S):'F12  5) 

221  F0RMAT(//1X'  BOMB  ATTITUDE  (AIRCRAFT  AXES)') 
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n 


i i 

23 

24 

25 
2b 

1 2 b 
2? 

28 

29 

30 

3 1 

32 

33 

37 

38 

39 

40 
42 

4 1 
C 

1 001 
1 002 


1003 

1004 

1005 
10  0b 
1 0 07 
1 008 
1009 
10  10 

C 

1 0 1 

c 

4002 

C 

C 


C 

C 

C 

1099 

C 

C 


20  1 


FORMAT! // IX ' ANGULAR  RATES  (BOMB  AXES)') 

FORMAT< IX'ROLL  RATE  < DE G /S  ) : ' F 1 2 . 5 ) 

FORMAT< 1 X ' P ITCH  RATE  < DEG/S > » ' F 12 . 5 ) 

FORMAT! 1 X ' Y AW  RATE  ( D EG / S ) : ' F 1 2 . 5 ) 

F0RMAT(//1X'***AIRCRAFT  PARAMETERS  •***') 

F0RHAT(/1X' HE  IGHT  ( HE TR E S ) : ' F 1 2 . 5 ) 

F0RMAT(/1X' ANGLE  OF  ATTACK  (DEG):'F12  5) 

FORMAT! /IX' DIVE  ANGLE  < D EG  ) : ' F 1 2 . 5 ) 

FORMAT! //lX'ANGULAR  RATES  (AIRCRAFT  AXES):') 
F0RHAT(//1X'***PR0BE  ATTITUDE  (TUNNEL  AXES  )*** ' ) 

FORMAT( IX'ROLL  ANGLE  ( 0 E G ) : ' F 1 2 . 5 ) 

FORHAT( lX'PITCH  ANGLE  ( D EG > : ' F 1 2 . 5 > 

FORHAT( IX' YAM  ANGLE  ( DE G > : ' F 1 2 . 5 ) 

FORMAT(//lX'***HISCELLANEOUS  INFORMATION***'  ) 

FORMAT( lX'CENTRIPETAL  ACCN  OF  AIRCRAFT  (IN  C ' S ) s " F 1 2 . 5 ) 
FORHATdX'TIME  STEP  (S):'F12.5) 

F OR M AT ( 1 X ' F RE E STREAM  MACH  N0.:'F12.5) 

FORMAT! /1X"B0NB-AL0NE  PRINT  INCREMENT  (IN  D T ' S ) : * I 4 ) 
FORHATdX'TIME  LIMIT  FOR  DROP  (S):'F12.5> 

F0RMAT(/1X'  ALTER  PARAMETERS?') 

FORMAT! / 1 X' ENTER  N. WHERE  N='/1X'1  B0MB;2  AIRCRAFT;3  PROBE; 
14  MISCELLANEOUS;' 

2 /1X'5  MEASUREMENT  STATIONS;b  DIAMETERS  AT  MS.'/) 

FORMAT! I 2 > 

F0RMAT(/1X' ENTER  NUMBER  OF  CHANCES') 

F0RMAT(/1X' XMS(L  ):  ENTER  L, CARRIAGE  RETURN. THEN  VALUE') 
F0RMAT(/1X' DMS(L  ):  ENTER  L. CARRIAGE  RETURN. THEN  VALUE') 
FORMAT(/lX' FILE  NAME  NOT  KNOWN  - TRY  AGAIN!'//) 

F0RMAT(/1X' ENTER  N , AS  DEFINED  ABOVE') 

F OR  M AT( IX':') 

F0RMAT(/1X' FILE  ALREADY  PRESENT  ON  DISK  - DO  YOU  WISH 
l TO  DISCARD  I T?  ( Y/N  )'  ) 

FORHATdX'CTYS  CONSTANTS  FOR  FILE  '.A5.A4//> 

FORMAT( A 1 ) 


TWOP 1=2 . *PI 
P I 4 = P I /4  . 

I OL  D = 1 
I PR  T D=  1 

ASK  FOR  FILE  NAME 

WRI TE( 4.1) 

READ(7.2)FILE 

ASK  IF  CREATING  NEW  FILE. 

I RS  T = 1 
CALL  CTRLP 

CO  TO  ( 201 , 202,203 , 204, 205 . 206, 207  ),  IRST 
WRI TE( 4, 3 ) 

READ(  7 ,4002  )ANS 

IF(ANS  EO  ANSNO)  GO  TO  1000 

CALL  FST AT( 2, F ILE,  I PRES  ) 

I F ( IPRES  EQ  0 > GO  TO  302 
WRITE! 4, 1010) 

READ( 7 , 40B2  ) ANS 
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I F( ANS . E Q . A NS  Y ES  > GO  TO  3B2 
CALL  C LO  S E(  2 ) 

GO  TO  1099 

C 

C YES.  NEW  FILE. 

302  I RS  T =2 

202  II R I TEC  4j  4 > 

C 

C ASK  FOR  POSITION  OF  6 OH  B C G.  (AIRCRAFT  AXES). 

URI TE<  4,  14) 

U R I T E ( 4/15)  X 
CALL  DECDE(X) 

UR  I T E(  4/16)  Y 
CALL  OECDE(Y) 

URITE(  4,  17)  Z 
CALL  DECDE(Z) 

ASK  FOR  VELOCITY  OF  BOHB  C.G.  (AIRCRAFT  AXES). 
URI  TE(  4,18) 

URI TE( 4,19)  U 
CALL  OECDE(U) 

URI  TE«.  4,  20  > V 
CALL  OECDE(Y) 

URI TE( 4,21)  U 
CALL  DECDE(U) 


ASK  FOR  BOHB  ATTITUDE  (AIRCRAFT  AXES). 
U R I T E(  4,121) 

U R I T E< 4, 31 >20LLS 
CALL  DECDE( ROLLS  ) 

UR  I T E( 4, 32)  PITCHS 
CALL  OECDE(PITCHS) 

URI TE( 4, 33)  YAUS 
CALL  DECDE(YAUS) 


ASK  FOR  ANGULAR  RATES  (BOHB  AXES). 

URITEt  4, 22) 

URI TE( 4,23)  PS 
CALL  DECDEl PS  ) 

URI T E( 4, 24)  QS 
CALL  DECDE( QS  ) 

URI  TE(  4,25)  RS 
CALL  DECDE( RS  > 

ASK  FOR  BOHB  PRINCIPAL  HOHENTS  OF  INERTIA 

URI  TE(  4,7  14) 

IX  ( KGH  H2  ) 

URI  TE(  4,5)  AX X 
CALL  DECDE(AXX) 

IY  (KGH  H 2 ) 

URI  TE(  4,7)  AY  Y 
CALL  DECDE(AYY) 

IZ  (KGH  H2 ) 

URI TE( 4, 8 ) A ZZ 
CALL  DECDE(AZZ) 

HASS  OF  BOHB  (KGH) 
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U R 1 TEC 4»  9 ) ST  HASS 
CALL  DECDE( STMASS > 

C 

C C.G.  POSITION  (METRES  FROM  TAIL) 

URI TE< 4, 10)  XCG 
CALL  DECDE(XCC) 

C 

C SCALE. AS  IN  1 /SCALE  . 

HRITE(  4.11)  SCALE 
CALL  DECDE( SCALE  > 

GO  TO  (3B3, 73  >, IOLD 

C 

C GET  THE  ’ NMS'  MEASUREMENT  STATIONS. 

3 B 3 I RS  T =3 

2B  3 HRI  TE(  4 , 1 2 )NMS 

DO  5BB  L = 1 . NMS 
SBB  CALL  DECDE<  XNS( L ) ) 

C 

C GET  DIAMETERS  AT  THESE  MEASURED  STATIONS. 

I  RS  T =4 

2 B 4 MRITE(4.  13)  NMS 

DO  5 B 1 L = 1 < NMS 
5 B 1 CALL  DECDE( DMS(L  )) 

C 

C*******ACQUIRE  AIRCRAFT  PARAMETERS****** 

C 

3 B 5 I RS  T = 5 

2 B 5 W R I T E<  4.  26) 

C 

C AIRCRAFT  HEIGHT  (METRES) 

URI  TE(  4,126)  HGHTB 
CALL  DECDE(HGHTB) 

C 

C ANCLE  OF  ATTACK  (DEGREES). 

HRI TE( 4, 27  ) ATTACS 
CALL  DECDEi ATTACS) 

ATTACK=ATTACS*DTR 

L 

C DIVE  ANGLE  (DEGREES). 

H R I T E( 4. 28)  DIVES 
CALL  DECDE( DIVES) 

DIVE=DIVES*DTR 

C 

GO  TO  ( 3B6, 73  ),  IOLD 

C 

C*******ASK  FOR  PROBE  ATTITUDE****** 

C 

3 B 6 I RS  T =6 

2 B 6 BRITE( 4, 3B) 

C PROBE  ROLL  (DEGREES) 

URITE(4.31)  SPTROL 
CALL  DECDE( SPTROL  ) 

APTROL=SPTROL*DTR 

c 

C PROBE  PITCH  (DECREES). 

HRITE( 4, 32)  SPTPIT 
CALL  DECDE( SPTPIT) 

APTPIT=SPTPIT*DTR 
CO  TO  ( 3B7. 73  ). IOLD 
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C 

c 

t,*****  acquire  miscellaneous  information***** 

c 

387  I RS  T = 7 

2 B 7 WRITEC  4 / 37) 

C 

C CENTRIPETAL  ACCN  OF  AIRCRAFT  (IN  G'S). 

MRI TE< 4, 38)  GF 
CALL  DECDECGF) 

C 

C T IHE  STEP  < SECONDS  ) 

WRITE ( 4i 39)  DT 
CALL  DECDE(DT) 

C 

C FREE  STREAM  MACH  NUMBER. 

WRITE ( 4<  48)  RHACH8 
CALL  DECDE<  RMACHB  ) 

NUMBER  OF  TIME  INCREMENTS  BEFuRE  PRINTING 
(DURING  BOMB-ALONE  TRAJECTORY). 

WRI TEC  4,42)  NPRINT 
CALL  I DECDi NPRINT) 

TIME  LIMIT  FOR  DROP  (SECONDS). 

WRI TE( 4, 4 1 > TL  IM  IT 
CALL  DECDEC TLIHIT) 

GO  TO  ( bb6b , 73  ),  IOLD 

C 

c ********* 

C ********  * 

C*«*******MODIFY  EXISTING  CONSTANTS  BLOCK**** 


1 0 50 

72 

73 

i 1 00 

1181 

3 B 9 


L.  : 


l OL  D =2 

CALL  F 3 T A T ( 2, FILE,  IPRES  ) 

I F ( IPRES .NE  B)  GO  TO  72 
WRI TEC  4,  1887) 

GO  TO  1 B 9 9 

CALL  3EEKC2.FILE) 

REA  D ( 2 ) T , IT  , TS 
CALL  CL0SEC2) 

CALL  CTRLP 
WRITEC  4 , 1081) 

READC  7 , 4002  > ANS 

IFCANS  Efl  ANSNO)  GO  TO  3 SB 6 

GO  TO  < 1 IBB, 1 181  ),  IPRTD 

I PR  T D=  2 

WRITEC  4,  1802) 

WRI TEC  4,  1 BB8  ) 

READC  7 , 1 B 83  ) NSYS 

IFC NSYS  LT. 1 . OR  NSYS. GT . 6)  GO  TO  1188 
GO  TO  (382,385,386,387, 3B9, 318), NSYS 
WRI TEC  4,  1 BB4  ) 

READC  7 , 1BB3  > NCH 

IFCNCH  LE  0 OR  NCH  GT.NHS)  GO  TO  73 
WRI TEC  4,  1 BBS  ) 

D 07  4 J J= 1 , NCH 
WRITEC  4,  1889  ) 

READC  7 , 1 BB3  ) L 
CALL  DECDEC XMSC L ) ) 

GO  TO  73 


' 
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C 

310  URITE( 4, 1004) 

READ( 7 . 1 003 ) NCH 

IF<NCH  LEO  OR  HCH.GT.HHS)  GO  TO  73 
UR  I TE<  4 < 1006) 

D 07  5 J J = 1 / NCH 
URITEC  4 . 1009  ) 
ft  E A D < 7 . 1003  ) L 
7 5 CALL  DECDE( DMS(L  )) 

GO  TO  73 
C 
C 

C FIND  EULER  ROLL/  PITCH  fc  YAU. 

C 

C 

6666  ROLL=ROLLS*DTR 

P ITCHR=PITCHS*DTR 
TPC  H = T AN( PITCHR) 

YAUSR=YAUS*DTR 
T YU  = TAN< Y AU  SR ) 

SRL=SIN< ROLL > 

CRL  = COS<  ROLL  > 

TANTH=TPCH*CRL+TYU*SRL 
P I T C H=  AT  A N( TANTH ) 

COSTH=COS(PITCH> 

TANB=<  TYU*CRL-TPCH*SRL  >*COSTH 
Y AU  = AT  AN< TANB  ) 

C 

C CALCULATE  FREE  STREAM  VELOCITY 

C 

TEHP=1 4 99-0  0 06 5* H GH T0 
VSOUND=SQRT( 40  1 742*TEMP+  109809 .48  > 
VEL I NF=VSOUND*RMACHB 


7 0 0 

C 

c 

c 


CALCULATE  AIRCRAFT  ANGULAR  RATES 
ASSUMES  A RO LL  = P I AND  AYAU*0 

P A=  0 

QA=GF*GRAVAC/VEL INF 
R A=  0 

CONVERT  TO  DEG/S  FOR  PRINTING. 


P AS  = 0 . 

OAS=OA*RTD 
R AS  = 0 

FIND  AREAS/ MAX  AREA. MAX  DIAMETER 

D MA X =0  . 

A MA  X = 0 

DO700  L= 1 , N MS 
A MS  < L > = P I 4*DHS( L >**2 
IF<  DNS(L  > LE  DMAX)  GO  TO  700 
0 MA  X *DHS  < L ) 

AMbX=AMS( L ) 

CONTINUE 


INITIALIZE  DIRECTION  COSINES  OF  B0N8  UR 


. 


T AIRCRAFT 


o o o o o 
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C 

SR=S IN( ROLL ) 

CR=COS(ROLL ) 

SP=SIN(PITCH> 

CP=COS( PITCH) 

sy=sin( yan > 

CY=COS( YAH) 

C 1 1 =CP*CY 

C 12=$R*SP*CY*CR*SY 
C 13=-CR*SP*CY*SR*SY 
C21 =-CP*SY 
C22=-SR*SP*SY+CR*CY 
C23=CR*SP*SY+SR*CY 
C'  3 1 = S P 
C 32= -SR* CP 
C33=CR*CP 
C 

C CALCULATE  ANGULAR  RATES  OF  BOMB  (IN  BOMB  AXES)  y R.T  FIXED 
C EARTH  AXIS  SYSTEM 

C 

C FIRST  CALCULATE  DELP,  DELfl  l DELR  FROM  PROJECTED 

C ROLL,  PITCH  AND  YAH  RATES. 

C 

RR=PS*DTR 

PP=fiS*DTR/COS( PITCHR)**2 
YY  = RS*DTR/COS<  YAUSR  )**2 

C 

DELP=C21  *(C32*YY-C33*PP  >+RR*CP/CY 
DELfi=Cll*(C33*PP-C32*YY  ) 

DELR  = C1 1*(C22*YY-C23*PP  ) 

NON  ADD  ( C ) . ( P A > » THE  AIRCRAFT'S  ANGULAR 
RATES  y R.T.  FIXED  AXES,  TRANSFORMED  TO  BOMB  AXES. 

< P>  = <C  ) ( P A > ♦ (DELP  ) 

D 08  0 1 1=  1 ,3 
S UH  = 0 . 

DO80B  J=l  , 3 
800  SUM  = SUM+C<  I ,J  >*P AV(  J ) 

80  1 PV<  I )=SUM+DELPV(  I ) 

C 

C 

C INITIAL  ANGULAR  MOMENTUM  OF  BOMB  (KGM  M2/S > 

C 

HX= AXX*P 
KY=AYY*fl 
HZ=AZZ*R 

C 

C CALCULATE  AIRCRAFT  ACCELERATIONS  (IN  AIRCRAFT  AXES). 

C 

CFG=GF*GRAVAC 
AX=GFG*SIN( ATTACK) 

AY=B 

AZ=-GFG*COS( ATTACK  ) 

C 

C NRITE  T BLOCK  ON  DISK 

C 


CALL  ENTER  (2, FILE) 
NRI TE( 2)  T, IT, TS 
CALL  CLOSE( 2 ) 
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PRINT  ALL  CONSTANTS  ON  THE  LINE  PRINTER 


MRI TEC  9, 6 74  9 ) T 

C 6 74  9 FORNAT< 1 X 10G1 3 5 ) 

C 

WRI TE< 9 , 101  ) FILE 
MRITEC9,  4 ) 

MRI TEC 9.  5 > AX X 
MRI TEC  9, 7 ) AVY 
MRI TEC  9, 8 > AZZ 
MRITEC9.9)  STNASS 
MRI TEC  9,  IB)  XCG 
MRI TEC  9/  1 l ) SCALE 
MRI TEC  9.14) 

MRI TEC  9.15)  X 
MRI TEC  9 . 18)  Y 
MRI TEC  9.  17)  Z 
MRI TEC  9.  18) 

MRI TEC  9.  19)  U 
MRI TEC  9. 2B)  V 
MRI TEC  9. 21  ) M 
MRI TEC  9.121) 

MRI TEC  9. 31  ) ROLLS 
MRITEC  9. 32)  PITCHS 
MRITEC9.33)  YAMS 
MRI TEC  9. 22) 

MRI  TEC  9.  23  > PS 
MRITEC  9. 24)  QS 
MRITEC  9, 25)  RS 
MRI TEC  9,  1B25 ) XHS 

1 025  F0RHATC71X' MEASUREMENT  STATIONS  A RE : ' / 1 X 1 2F 1 0 . 4// > 

MRITEC  9.  1026)  DNS 

1028  FORNATC IX' DIAMETERS  AT  THESE  STATIONS  ARE : ' / IX 1 2F 10 . 4/V ) 
MRI  TEC  9,  26  ) 

MRITEC  9,126)  HGHT0 
MRITEC 9, 27)  ATTACS 
MRITEC 9, 28)  DIVES 
MRI TEC  9, 29  > 

MRI TEC  9, 23)  PAS 
MRI TEC  9, 24  ) Q A S 
MRI TEC  9. 25  > RAS 
MRI TEC  9, 3B  ) 

MRITEC 9. 31)  SPTROL 
MRITEC  9, 32)  SPTPIT 
MRI TEC  9, 37  ) 

MRITEC  9, 38)  GF 
MRI TE< 9.39)  DT 
MRITEC  9, 40  )RHACHB 
MRI TEC  9,42)  NPRINT 
MRI TEC  9,41)  TL  INIT 
STOP 
END 
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C DECDE  DATE:  12-7-77 

C 

C ROUTINES  TO  ACQUIRE  A FLOATING  POINT  OR  FIXED 
C POINT  NUMBER  FROM  THE  TELETYPE. 

C IF  THE  OPERATOR  HAS  TYPED  A CARRIAGE  RETURN  ONLY , 

C THE  VALUE  OF  THE  ARGUMENT  IS  NOT  CHANCED  AND  AN  IMMEDIATE 
C RETURN  TO  THE  CALLING  PROGRAM  IS  EFFECTED. 

C 

C CALLED  BY  EDCON  . 

C 

C 


SUBROUTINE  DECDE  (A) 

DIMENSION  ASC( 3 ) 

DATA  BLAHKS/5H  / 

C 

1*1 

GO  TO  IBB 

C 

ENTRY  IDECD  <N> 

1=2 

IBB  READ<  7 , 1 > ASC 

IF( ASC< 1 > EQ. BLANKS > GO  TO  3BB 
GO  TO  < 1B1, 1B2  >,  1 
1 B 1 DECQOEC 12. ASC  , > A 

CO  TO  3BB 

1 02  DECODE( 12, ASC,  > N 

1 FORM AT( 3 A5 ) 

3 B B RETURN 

END 


i 


1 
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APPENDIX  II 

LISTING  OF  THE  MAIN  PROGRAM 


C HhIN  DATE:  8-7-77 

C 

••♦•♦PROGRAM  MAIN  (MIND  TUNNEL  VERSION) 

•••••ROUTI NES  REQUIRED  . 

TRANS, PAR A8.HRCT AN- CPBIN, READ  IN, TRFREE. TRSET 
FLOW, CALC, TRVCTS, PROBE 
CREF , C INC , CORECT 
C I NT , DAUX 

C ***** L 1ST  OF  STOP  LOCATIONS... 

C * ** •*  l -PROBE, 2-READIN,3-, 4-, 5-PROBE, 6-READIN 
C* ****? -REA  DIN, 8-, 9-,  IB- CREF,  11- PROBE. 

C 

DIMENSION  F ILEC 2 ), TFILE< 2 ) 

DIMENSION  TMACHF(2  ) ,PRATFL<2  ),FILEMC(2  ),ST0RCX<2  > 

DIMENSION  A MK  8 2<  2 ) , AK  AR(  2 ) 

DIMENSION  F YEC<  6 ),  F SV ( 3 , 6>,TIMES(  3 > » AB  C ( 3 , 8 ) 

DIMENSION  XF(3),XN(3),DUM(3),ELA<3) 

C 

C 

C***** STORAGE  FOR  AREAS, DIAMETERS, DYNAMIC  PRESSURE, MACH  NUMBER, 
C***** VELOCITY  U , V , U , MEASUREMENT  STATIONS  (AT  EACH  MS, BOMB  AXES) 
DIMENSION  AMS<  12  ),  DMS<  1 2 ),  QMS<  12),  RMS<  12), 

1 UMS<  12  ),  VMSC  12  ),  UMS<  1 2 ),  XMSC  12) 

C 

C 

COMMON  T ( 320),  IT(5  >,NF 

C 

C 

C**.**TIME, INTEGRATION  TIME  STEP 

EQUIVALENCE  ( T < 2 ) , T I M E ) , ( T < 3 ) , D T ) (AIRCRAFT 

C*»***p  IRECT  ION  COSINES  DEFINING  ORIENTATION  OF  BOMB  HITH  RESPECT  TO 
EQUIVALENCE  ( T < 4 ) , C 1 1 ) , < T<  5 ) , C2 1 ) , ( T< 8 ) , C3 1 ) 

EQUIVALENCE  ( T<  7),C12),<T<  8),C22),(T<  9),C32) 

EQUIVALENCE  < T ( 1 B ) , C 1 3 ) , < T ( 1 1 ) , C2 3 ) , < T < 1 2 ) , C 33  ) 

C*****PQS IT  I ON  OF  BOMB  CENTRE  OF  GRAVITY  IN  AIRCRAFT  AXES 
EQUIVALENCE  <T<13),X),<T<14).Y),(T<15),Z> 

C * **  **  V EL  OC I TY  OF  BOMB  C.G.  IN  AIRCRAFT  AXES. 

EQUIVALENCE  <T<16).U),<T<17).V),<T<18).M) 

C*****ANGULAR  MOMENTUM  OF  BOMB  IN  BOMB  AXES 

EQUIVALENCE  < T<  1 9 ) , HX  ),  ( T<  20  ) , HY  ) , < T(  2 1 ),  H2  ) 

C*****ratE  OF  CHANGE  OF  DIRECTION  COSINES  ABOVE 

EQUIVALENCE  ( T ( 2 2 > , DC  11  ) , ( T<  2 3 ) , D C 2 1 ) , < T( 24  ) , DC3 1 ) 
EQUIVALENCE  < T( 25  ) , DC  1 2 ) . ( T<  26  ) , DC22  ), ( T<  27  ) , DC32  ) 
EQUIVALENCE  < T ( 2 8 ) , DC  1 3 ) , < T( 29  ) , D C 23  ) , < T( 3B  ) , DC3 3 ) 

C*****RATE  OF  CHANGE  OF  X,Y.Z, ABOVE 


EQUIVALENCE  ( T(  3 1 ) , DX  ),  ( T<  32  ) , DY  ) , C T<  33  ).  DZ  ) 
C*****RATE  OF  CHANGE  OF  U.V,H  ABOVE 

EQUIVALENCE  < T<  3 4 ) , DU  ),  < T<  35  ) . DV  ) , < T<  36  ),  DU  ) 
C*****RATE  OF  CHANCE  OF  HX,HY,HZ  ABOVE 

EQUIVALENCE  < T < 3 7 ) , DH X ) , ( T C 3 8 ) , DH Y > . ( T < 39  > , D HZ  ) 


C*****T(40)  TO  T ( 75  ) RESERVED  FOR  INTH. 

C***** A IRCRAFT  ROTATION  ANGLES  M.R.T  EARTH  AXES 

EQUIVALENCE  ( T < 7 6 ) . AR OL L ) , < T ( 77  ) . A P I TC H ) . ( T ( 78  ) , A YA M ) 
C*****80MB  MOMENTS  OF  INERTIA 

EQUIVALENCE  ( T(  79  ) , AXX  ) . < T(  80  ),  AYY  ),  ( T<  81  ),  AZZ  ) 
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C*****BOHB  NASS  ( KCM  SLUGS) 

EQUIVALENCE  < T( 82 ) , STNASS  ) 

C * «***P  ROBE  ROLL  AND  PITCH  ANGLES  M.R.T.  TRAVERSE  AXES 
EQUIVALENCE  < T( 83  ) , APTROL  > , < T< 84  ) . APTP I T ) 

C**»**T  IHE  LIMIT  FOR  DROP, AIR  TEMPERATURE  (DEG  C>, 

C *****  B OM  B C.G.  (METRES  FROM  TAIL), AIRCRAFT  HE  I GHT( METRES >, 
C*****BOMB  SCALE  (AS  IN  1/SCALE  ). T( 98  ) SP ARE , SPEED  OF  SOUND/ 
C*****ACCELERATION  DUE  TO  GRAVI TY. T( 93  ) SPARE 

EQUIVALENCE  ( T ( 85 ) . TL I M I T ) , ( T ( 86 > . TEMP  ) , ( T( 87 > , X C G > / 

1 ( T(88>, HGHTB  ),( T( 89  ), SCALE >, 

2 ( T(91 >, VSOUND),( T(92>, GRAVAC) 

C *****  TOT  AL  PITCH  AND  ROLL  OF  BOMB  IN  FLOHFI ELD  AT  REFERENCE  POINT. 

EQUIVALENCE  ( T ( 9 4 ) , THET OT  ) , ( T( 95  ) , PH  I TOT  ) 

C ***** F RE E STREAM  MACH  NUMBER, MAX  BOMB  AREA, MAX  BOMB  DIAMETER 
EQUIVALENCE  ( T ( 9 6 > , RM AC H B ) , ( T ( 97  ) , AM AX  ) , ( T( 98), DM AX) 
c*****aircraft  ANGLE  OF  ATTACK 

EQUIVALENCE  ( T( 99 ), ATTACK  ) 

C * **  **  T ( 1 BB  > TO  T(1B3)  SPARE. 

C 

C*****FR£E  STREAM  DYNAMIC  PRESSURE 
EQUIVALENCE  ( T( 1 84  ) , QDPB  ) 

C *****  u IND  FORCES  ON  BOMB  IN  AIRCRAFT  AXES. 

EQUIVALENCE  ( T ( 1 B5  ) , UX > , ( T ( l B6 ) , U Y ) . ( T( 1 B7  ) . UZ  ) 

C *****  V IND  MOMENTS  ON  BOMB  IN  B0N8  AXES. 

EQUIVALENCE  (T(1H8>,«L),(T(1B9),WM).(T(I1B>.HN> 

C * **** A IRCR AFT  R OLL , P I TC H , Y AN  RATE  ABOUT  AIRCRAFT  AXES 
EQUIVALENCE  ( T ( 1 1 1 > , P A ) , ( T ( 1 1 2 ) , Q A ) , ( T ( 1 1 3 ) . RA  ) 

C * ****  FREE  = -1  (M.)  IF  BOMB-ALONE  TRAJECTORY 

C IS  NOT  (IS)  OPERATING  . 

EQUIVALENCE  ( T ( 1 1 4 > . F RE E > 

C * ****ALPHA  AND  BETA  OF  BOMB  USED  FOR  CALCULATING  COEFFIC  ENTS 
EQUIVALENCE  ( T ( 1 1 5 ) , ALPH A > , ( T( 1 1 6 ) , BET A > 

C***** AIRCRAFT  ACCELERATIONS  IN  AIRCRAFT  AXES. 

EQUIVALENCE  ( T( t 1 7 ) , AX  ) , ( T( 1 1 8 > , A Y ), ( T(  1 1 9 > . A2  ) 

C *****  F OR  CE  S ON  80MB  C.G.  IN  AIRCRAFT  AXES. 

EQUIVALENCE  ( T( 1 2B  ) , FX, F VEC ) , ( T( 1 2 1 > , F Y ), ( T( 1 22 > , FZ > 

C***** MOMENTS  ON  BOMB  ABOUT  BOMB  X,Y,Z  AXES 

EQUIVALENCE  (T(123),FL),(T(124),FM),(T(125),FN) 

C*****ROTAT I ON  ANGLES  DEFINING  ORIENTATION  OF  BOMB  U . R . T A I RC R AF T . 

EQUIVALENCE  ( T ( 1 26  ) , R OL L ) , ( T ( 1 2 7 ) , P I TC H ) . ( T ( 1 2 8 ) , YAM) 
C*****ROLL, P ITCH . YAN  RATES  OF  BOMB  ABOUT  BOMB  X.Y,2  AXES 
EQUIVALENCE  ( T ( 1 29  ) , P ),  ( T(  1 3 B ) . Q ) , ( T(  1 3 1 ) , R ) 

C*****pRESSURE  AT  HEIGHT, DYNAMIC  PRESSURE  AT  HEIGHT. 

EQUIVALENCE  ( T( 1 32  > , PRE SS  ) , ( T( 1 33  > , QDP ALT  ) 

C*****DECREES  TO  RADIANS, RADIANS  TO  DEGREES 

EQUIVALENCE  ( T( 1 34  > , DTR  ) , ( T( 1 35  ), RTD  ) 

C * **** VELOC I TY  OF  AIR  AT  -INFINITY- 
EQUIVALENCE  ( T( 136 > , VELINF  ) 

C*  **** A IRCRAFT  SELF-ROTATION-RATE  MATRIX. 

EQUIVALENCE  ( T ( 1 37 ) , A R 1 1 ) , ( T ( 1 3 8 ) , AR 2 1 ) , ( T( 1 39  > , AR3 l ) 
EQUIVALENCE  ( T( 1 4B  ) , AR t 2 ) . ( T( 1 4 1 ) , AR22  ) , ( T( 1 42  ) , AR32 ) 
EQUIVALENCE  ( T ( 1 43  ) , AR 1 3 ) , ( T ( 1 44 ) , AR23  ) , ( T( 1 45  ) , AR33 > 

C * **  **  8 0M8  SELF-ROTATION-RATE  MATRIX. 

EQUIVALENCE  ( T( 1 46  ) , BR 1 1 ) , ( T( 14 7 ) , BR21  ) . ( T(  1 48  ), BR3 1 ) 
EQUIVALENCE  ( T ( 1 49  ) , B R 1 2 ) , ( T ( 1 5B  ) , BR 22  ) . ( T( 1 5 1 ) , BR3 2 ) 
EQUIVALENCE  ( T( 1 52  ) . B R 1 3 ) , ( T ( 1 5 3 ) , BR23  ) . ( T( 1 54  ) , BR3 3 ) 

C* **** T( 1 55  ) TO  T( 1 6B  ) SPARE. 


WSR1.-0005-TR 


48 


C * **  **  G RA  V I T Y FORCE  IN  AIRCRAFT  AXES. 

EQUIVALENCE  <T<i6l),GX).<T<162).CY>»<T<163>,CZ) 

C*****ORI GI NAL  DIVE  ANGLE  OF  AIRCRAFT. 

EQUIVALENCE  < T< 164 >.DIVEB) 

£*****PI, TyoPI /CONVERSION  FROM  METRES  TO  INCHES 

EQUIVALENCE  ( T ( 1 65 ) , P I ) . < T< I 66 > . THOP I > - ( T< I 67  ) . C H I ) 

C CENTRIPETAL  ACCELERATION  OF  AIRCRAFT  (IN  G‘S>. 

EQUIVALENCE  ( T< 168  >,GF> 

C*****AT  ALL  MEASUREMENT  STATIONS  . . AREAS, DIAMETERS, DYNAMIC  PRESSURE. 
C*****MACH  NUMBER  , U , V , Q (WIND  VELOCITY  COMPONENTS . BOMB  AXES). 

C * **  **  A ND  MEASUREMENT  STATIONS  <C  G.*B,NOSE  +VE  . I . E . . BOMB  AXES) 
EQUIVALENCE  ( T< 1 69  ) . AMS  ) , ( T( 1 8 i ) , DMS  ) . ( T< 1 93  ) . QMS  ) , 

1 < T<205  ) ,RMS  ).<  T<  217  ).  UMS  ),  < T<  22  9 ).  VMS  ),  < T<  24  i ).  QMS). 

2 < T(253  >,XMS  ) 

C*****FREE  STREAM  AERODYNAMIC  COEFFS  BOMB  AXES 

EQUIVALENCE  ( T< 265  > . CX  ) . < T< 266  ) . C Y >. ( T<  26  7 > . CZ  ). 

1 (T(268)/CL)/(T<269).CM)i<T<2?B).CN) 

C*****  INCREMENTS  TO  COEFFS  DUE  TO  NON-UNIFORM  FLOW 

EQUIVALENCE  < T < 2 7 1 > . D CX  ) . ( T(  2 72  ) . D C Y ) . ( T(  27 3 ) . DC Z ) . 

1 < T<  2 74  ) ,DCL  ) , ( T<  2 75  ),  DCM  ).  < T<  27  6 >.  DCN  ) 

C *****  pro  BE  PRESSURE  MEASUREMENTS 

EQUIVALENCE  < T ( 2 77  > , D P 1 3 ) » < T ( 27 8 ) . DP  24  ) . 
l ( T<  2 79  > , PHAN  ). < T<  28 B) , PPIT  ).<  T<  281  ). T PRESS  ) 

C AREA  TO  SAVE  COEFFICIENTS  AT  PREVIOUS  2 TRAJ  . POINTS 

EQUIVALENCE  ( T(282).FSV  ).(T<  3BB). TIMES  >.(T(3B3).ABC  ) 

C 

C 

C *****  N 0 OF  MEASUREMENT  STATIONS, REF  PO I N T , I NT E GR AT  I ON  LIMITS. 
EQUIVALENCE  ( I T(  1 ) , NH S ) . < I T<  2 > , NR E F ) , ( I T<  3 ) . NB AC X ) 

EQUIVALENCE  < I T< 4 ) , NFRONT  ) , < I T< 3 ) , NPRI NT  ) 

C 

EQUIVALENCE  < XF( 2 ) , YF  >, ( XF( 3 > , ZF  ) 

EQUIVALENCE  ( X N<  1 ) , XN 1 > . < X N<  2 ) , YN  ) , < XN < 3 ) , Z N ) 

C 


DATA 

C 

C*»***L ABELS 

C 

DATA 

DATA 

DATA 

C 

C***** LAB  ELS 

C 


C 

C 


DATA 

DATA 

DATA 


TFILE/5HTFILE, 4H2CTS/, ANSYES/SHY 

FOR  YAWMETER  PROBE  FILES 

TNACHF/5HTMACH ,4HTPRB/ 
PRATFL/5HPRATI .4H0PRB/ 

F ILEMC/5HMACHC ,4HXPR8/ 

FOR  STORE  COEFFICIENT  FILES 

AMK82/'5HHK82C,  4HXSTC/ 
AKAR/5HKARCX, 4HSSTC/ 

ANSM/SHH  / , ANSK/5HK  / 


C 


T AN  ( X ) = S I N< X >/COS( X > 


I MP  A C T = 1 

C 

C SET  UP  FOR  INTERRUPTS  FROM  TRAVERSE  RIG 

C READ  8UT  TON 

C 


/ 


CALL  TRSET  ( I FREE  > 
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C 

C 

C SET  UP  THE  PROBE  PILES  FOR  DIRECT  ACCESS  I/O 

C 

CALL  DEFINE! 14,8, 1958. TNACHF.  IVl.B.B.fl) 

C HAS  28 15  = 85X3 1 ; NOW  65X3B  <PITCH=3B  ABSENT) 

C 

CALL  DEFINE! 3, 8 , 33B , PRATFL , I V2, B, 8 , B ) 

C WAS  341=31X11;  NOW  3BX1 1 (PITCH-3B  ABSENT) 

C 

CALL  DEFINE! 12, IB. 18491 , FILEMC, I V3,B,B, B> 

C 

C 

WRITE! 4, 6783) 

6783  F0RHAT(///1X’ CAPTIVE  TRAJECTORY  YAWNETER  SYSTEM'/) 
WRITE! 4, 6784) 

6784  FORH AT! / 1 X' HEW  TRAJ?  !Y  OR  N)') 

READ<7,6?95>  ANS 

6795  FORHAT(Al) 

IF! ANS  £Q  ANSYES  ) CO  TO  8888 
CALL  SEEK!2,TFILE> 

READ12  )T , IT ,STORCX,  INDST ,HF 
CALL  CLOSE! 2 ) 

GO  TO  8496 

8888  WRI TE! 4, 66B1  ) 

6681  FORMAT! /IX' Hi!  82  !H)  OR  KARINGA  !K)?') 

READ!7.6795>  ANS 

IF  ! ANS . HE . ANSH)  CO  TO  6602 

STORCX! 1 >»AHK82!  1 > 

STORCX! 2 )»AHK82! 2) 

CO  TO  8889 

66B2  IF  !ANS  HE  ANSK)  GO  TO  8888 
STORCX! 1 )=AKAR( 1 ) 

ST0RCX(2  >*AKAR<2  > 

8889  WRI TE!  4, 6796  ) 

6796  F0RHAT!/1X'ENTER  FILE  NAME') 

6782  READ!  7,6794  ) FILE 

6794  FORMAT! A5.A4) 

CALL  FSTAT! 2, FILE,  IPRES  ) 

IF!  IPRES  HE  8 > GO  TO  67B6 
WRITE! 4, 6785) 

6785  F0RMAT!/1X' FILE  NOT  PRESENT  - TRY  AGAIN!'  ) 

GO  TO  6782 

6786  CALL  SEEK!2,FILE) 

READ! 2 )T , IT 

CALL  CLOSE! 2 ) 

I HD  S T=  1 
FRE  E *- 1 
C 

8496  CALL  DAUX 

GFG  = GF  *GR AVAC 
C *****  S TART  THE  FLIGHT* 

C 

3 DIVE=DIVEB-CFG*TIHE/VELINF 

C 

C AIRCRAFT  OR IEHTATI OH  W.R.T.  EARTH  AXES. 
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C 

A R0  L L*  P I 

APITCH«-(DIVE-ATTACK> 

C 

C NOTE:  ftPITCH  IS  ftN  EULER  ANGLE/  WHERE  WE  NAVE 

C ROLLED  THROUGH  PI.  HENCE.  THE  ABOVE 

C FORMULA  IS  CORRECT  EVEN  THOUGH  IT  LOOKS 

C AS  THOUGH  IT  SHOULD  HAVE  A MINUS  THROUGH  IT! 

C 

AYAW»B . 

C 

IF  < FREE  GT  B ) GO  TO  4 
IF  (ABS(CF).GT.l  . £ - 6 ) CO  TO  941 
SV3 VEL INF*T  IME 
VGT=B . 

CO  TO  942 

941  VC=2  *VELINF/GFG 
VGT  = T I ME / VC 
SV=VG*VELINF*SIN(VGT) 

942  HEIGHT»HCHTB-SV*SIN(DIVEB-VGT  > 

I F ( IFREE . EQ  B > GO  TO  5 

CALL  T RF  REE  < S TO  RC  X , H E I C HT  > 

I HD  S T 3 1 

I PR  I NT  =NPRI NT  - 1 

4 HEI GHT -Z 

C**»**CALCULATE  QUANTITIES  DEPENDING  ON  HEIGHT  OF  AIRCRAFT  (METRES). 
C***»*STATIC  PRESSURE  (NEWTONS/SQ  H> 

5 PRE  S S= 1 B 3 35  1 1 *GRA VAC*( 1 . -B . BBBB22359*HEI GHT  >**5 . 23S1 B3 
C*****TEMP£RATURE  (DEG  C) 

TEHP«14  99-B . BB63*HEIGHT 
C * **  **  V EL  OC  I TY  OF  SOUND  (METRES/SEC) 

VSOUND3SQRT ( 4B1 . 742*TEMP*1B98B9 .48  ) 

C 

IF  (FREE  GT  8 > RNACHB*SQRT(U**2*V**2  + W**2)/'VS0UHD 
C *****  v EL  OC  ITY  AT  INFINITY  (M/S) 

VEL INF=VSOUHD«RMACHB 
C*****DYNAMIC  PRESSURE  (HEWTOHS/SQ  M) 

QDPALT=B  5*1  4*PRESS*RMACHB**2 
QAALT=QDPALT*AMAX 
Q AD  3 QA AL  T *D  M A X 

C *****  F IND  AIRCRAFT  QUANTITIES  IN  AIRCRAFT  AXES. 

C 

C 

D UM ( 1 )>B  . 

D UM < 2 ) 30 
C 

IF  ( FREE  GT  B . > GO  TO  IB 

C 

C***.*FIND  GRAVITY  IN  AIRCRAFT  AXES. 

C***«*CRAVI TY  FORCE  ON  BOMB  IN  EARTH  AXES  -(B.B.GZE) 
GZE*-STNASS*GRAVAC 
D UM ( 3)»CZE 

C *****  F IND  GRAVITY  FORCE  ON  BOMB  IN  AIRCRAFT  AXES  FRONl  EARTH  AXES 
CALL  TRANS  ( C X , AROL L , AP I TCH. A YA W, ♦ 1 , DUN  ) 

GO  TO  11 


o n n o o o o o o o oonr>o*-r>r>  r>  o o •-*  o o r> 


- 51  - 


WSRL-0005-TR 


•••••FIND  WIND  FORCES  IN  AIRCRAFT /MOMENTS  IN  BOMB  AXES. 

B IPRINT*IPRINT+1 

CONVERT  VELOCITIES  FRON  EARTH  AXES  TO  BOMB  AXES. 
B0H8-AL0NE  TRAJECTORY 
CALL  CTRANS  (DUH,*1,U> 

UHS<  HREF  )=-DUH(  1 ) 

V NS  ( HREF  >»-DUH( 2 ) 

WNS ( NREF  >*-DUH( 3 > 

VEL  REF  =VEL I NF 

RHSREF  =RH ACHI 

RNS<  NREF  >*RHACHB 

VBE  T A = VHS( NREF  >/UHS(NREF  > 

WALPHA  = WNS<  NREF  >/UHS( NREF  ) 
THETGT=ATAN(SQRT(VBETA**2*WALPHA**2)) 

PHI  TOT  =>ARCTAN( VBETA/W ALPHA  ) 

RAT  I 0* 1 . 

CALL  CREF 
WX8*CX*QAALT 
WY6=CY*QAALT 
WZB • CZ  *Q  A AL  T 
WL=CL*QAD 
WH=CH*QAD 
WN»CN*QAD 
GO  TO  12 

*****F1RST  FIND  FLOW  PROPERTIES  ALONG  BONB  IN  BOHB  AXES 

1 CALL  FLOW 


CALCULATE  EFFECTIVE  INCIDENCE  (ALPHA)  AND 
EFFECTIVE  SIDESLIP  (BETA)  FRON  FIN  CAHBER 
AND  REFERENCE  POINT  INCIDENCE  AND  SIDESLIP. 

CALL  CORECT(BETAiVNS) 

CALL  CORECTt ALPHA, WHS > 


•••••FIND  TOTAL  PITCH  AND  TOTAL  ROLL  OF 
REFERENCE  POINT. 

•••••(TOTAL  ROLL  (PHITOT)  IS  ANCLE  THAT 
•••••IS  ROTATED  FRON  TOTAL  AXIS  *Z) 
VB£TA=TAN(BETA> 

WALPHA»T AN( ALPHA  ) 


•••••FIND  TOTAL  PITCH  ANGLE 

THET0T=ATAN(SQRT(VBETA**2*WALPHA* 
•••••FIND  ROLL  ANGLE 

PHI TOT*ARCT AM( VBETA.W ALPHA  ) 
•••••FIND  FORCE  AND  NONENT  COEFFS  DUE  TO 
•••••TO  FLOW  AT  REFERENCE  POINT  IN  BONB 
VELREF=SORT(UNS( HREF >**2+VNS( HREF 
RNSREF  = RNS( HREF  ) 

CALL  DEFINE( 13,29, 1664,ST0RCX,IV4 
CALL  CREF 
CALL  CLOSE(  13  ) 

C*«***ADJUST  TABULATED  COEFFS  TO  CALCULAT 
C*****SANE  DYNAHIC  PRESSURE 
RAT  I 0=QNS( NREF  )/QDPB 


BONB  W R T TOTAL  AXES  AT 
BOHB  +Z  AXIS 


• 2)  ) 


UNIFORN  FLOW  EQUAL 
AXES 

>**2+WHS(NREF)«*2) 

, B,  B,  I ) 


ED  COEFFS  BY  BASING  ON 


jA 
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C 

C*****FIND  INCREMENTS  TO  ABOVE  DUE  TO  FLOW  VARIATIONS  ALONG  BOMB 
CALL  C INC 

C * **  **F  I ND  FORCES  AND  MOMENTS  FROM  COEFFS  AND  Q AT  ALTITUDE  OF  FLIGHT 
QAD*QAALT*DMAX 
UXB=CX*RATIO*OAALT 
WYB=(CY*RATI0*DCY>*8AALT 
UZB*(CZ*RATIO+DCZ)*QAALT 
UL=CL*RAT  I0*8AD 
MM=( CM*RATIG+DCM  )*6AD 
WH=  < CN*R ATI O+DCN  )*  8 AD 

C *****  E XP  RE  S S FORCES  IN  AIRCRAFT  AXES  FROM  B0H8  AXES 
12  CALL  CTRANS  (UX.-l.MXB) 

C****« UPDATE  BOMB  ORIENTATION  IN  PROJECTED  YAM, PITCH  PLANE  OF  AIRCRAFT 
ROLL=ARCT AN  ( -C32, C33  ) 

PITCH* ATAN< -C 13/C1 1 ) 

Y AW  = AT  AN< C 1 2/ C 1 1 ) 

C 

C FIND  PROJECTED  ROLL,  PITCH  fc  YAH  RATES. 

C 

R0LLD=(DC33*C3  2-DC32*C33  >/< C 32**2+C33**2 > 

PITCHD=<  DC  1 1*C13 -DC 13*C1 1 >/<  Cl  1**2 +C1 3**2  > 

YAMD*  ( DC12*Cll-DCll*C12>/<  C 1 1 **2*C 12**2  ) 

C 

C*****ADD  THE  FORCES  AND  MOMENTS  TOGETHER* 

FX=WXfGX 
FY=y Y+GY 
FZ=HZ*G Z 
FL=HL 
FH=HM 
FN*UN 
C 

TIMES( 3>»TIMES<2  ) 

TIMES( 2)  = TIMES( 1 ) 

T IMES<  1 > = T I HE 
D Ob  8 L * 1 , 6 
C 

C SAVE  COEFFICIENTS  AT  PREVIOUS  2 POINTS. 

C 

FSV( 3, L >*FSV< 2 , L ) 

F S V ( 2,  L = F S V < 1 ,L  ) 

F S V ( 1 , L ) = FVEC<  L ) 

I F< INDST  GT  2)  GO  TO  6? 

ABC ( 1 , L > = FVEC( L ) 

A BC ( 2. L>*0. 

ABC ( 3. L >*B . 

GO  TO  68 

b 7 CALL  PARAB  ( A B C< 1 , L > , T I M ES , F S V( 1 , L > ) 

6 8 CONTINUE 

I NDST* 1NDST* l 
CALL  DAUX 
ROLDEG«ROLL*RTD 
PITDEG*PITCH*RTD 
Y AMDEG*Y A W*RTD 
PDEG*ROLLD*RTD 
80EG=PITCHD*RTD 
RDEG*YAUD*RTD 


p 
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C 

IF  ( FREE  GT  . 0 . AND  IPRINT . ME . MPRIMT)  GO  TO  5431 
IPRINT-0 

C 

5 4 30  IF  ( FREE . GT  . 0 . ) GO  TO  69 

C 

C CONVERT  U < V i W INTO  PHYSICAL  VELOCITIES  OF  BOMB 
C C C.  RELATIVE  TO  AIRCRAFT  AXES,  BY  ADDING  (AROT).(X> 

C 

UPR=U*RA*Y-QA*Z 
VPR*V+PA*Z-RA*X 
yPR=M+9A*X-PA*Y 
GO  TO  5432 

C 

69  UPR  *U 

V PR  = V 
upR-y 

c 

5432  URITE( 9, 6901 > T I HE , X , Y , Z , R OL DEC , P I TDEG , YAMD E G , 

1 UPR, VPR/MPR, PDEG, QDEG/RDEG 
BRITEC  9, 6 97 3 > F X , FY , FZ,FL,FH, FN 
TH£TOD=THETOT*RTD 
PHI TOD=PH  ITOT*RTD 

NRI TE( 9, 687 3)  THE  TOD, PHI TOD,VELREF, RHSREF, RATIO 

C 

GO  TO  (5431,6907),  IMPACT 

C 

68  73  F OR HAT(5X'T HE  TOT, PHITOT,VELREF, RHSREF, RATIO!  ' , 2G  15 . 5 , 3G  1 5 . 6 

1 ///) 

6973  F0RHAT(5X'FX, FY,FZ,FL,FH,FN!  ',6015.6) 

6901  FORHAT(IX'TIHE,X,Y,Z',F8.4,3X,3G14.5,10X'ROLL,PITCH,YAM', 

1 3G145/6X'U,V,y,,llX,3G14.5,4X'RATES:R0LL,PITCH,YAM', 

2 3G  14  5/  ) 

C *****  I NT  EG  RATE  ONE  T I HE  STEP* 

5431  CALL  I NT  M 

CALL  ENTER( 2, TFILE > 

H RI T E<  2 ) T, IT , STORCX,  INDST,NF 
CALL  CLOSE( 2 ) 

C 

IF  (FREE  GT.0  AND  Z.LT.0  ) GO  TO  16 

C 

C ***** S EE  IF  REACHED  TIHE  LIMIT  FOR  DROP* 

IF  ( TIME . LT  . TL  IHIT  ) GO  TO  3 
C*****F  INISHED  FLIGHT 

C 

y RI TE( 9, 6905) 

69B5  F0RHAT(/1X' TIME  LIMIT  EXCEEDED') 

CO  TO  6907 

16  yRITE(9,6904) 

6904  F0RHAT(/1X' STORE  HAS  STRUCK  GROUND^/') 

IMPACT -2 
CO  TO  5430 

C 

6907  CALL  CLOSE( 3 > 

CALL  CLOSE(  12  ) 

CALL  CLOSE(  13  > 

CALL  CLOSE(  14  ) 

STOP 

END 
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APPENDIX  III 

LISTINGS  OF  TRSET  AND  TRFREE 


/TRSET  D ATE : 16-6-77 

B 

/ROUTINE  TO  SET  UP  THE  SKIP  CHAIN  FOR  INTERRUPTS 
/FROM  THE  TRAVERSE  RIG  READ  BUTTON,  AND  ALSO  TO 

/SERVICE  THAT  INTERRUPT 

$ 

i 

TRSF»7B500l 

/ 

GLOBL  TRSET,  DA 

/ 

/CALLING  SEOUENCE:  CALL  TRSET  (IFREE) 

/ 

/ 

TRSET  XX 

J NS  * .DA 
J HP  . + 2 
IFREE  0 

CAL 

16 

TRSF 
TRBINT 
DZH*  IFREE 
J HP  * TRSET 


/INTERRUPT  SERVICE  ROUTINE 

B 

TRBINT  DAC  ACSVE# 

LAC*  <0 
DAC  PCSVEI 

705002  /CLEAR  THE  FLAG. 

I SZ*  IFREE 
LAC  ACSVE 
ION 

J HP  * PCS VE 
END 


o o 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


TRFREE  DATE:  12-7-77 

ROUTINE  TO  CHANGE  F R OH  AN  AIRCRAFT-ORIENTED  TO  AN 
EARTH-AXIS-ORIENTED  TRAJECTORY. 

OPERATES  AT  THE  POINT  WHERE  THE  STORE  IS  CONSIDERED  TO  HOVE 
FROH  A TRAJECTORY  WHICH  IS  INFLUENCED  BY  THE  PRESENCE  OF 
THE  AIRCRAFT  TO  A UNIFORH-FLOW  TRAJECTORY. 

SUBROUTINE  TRFREE  ( ST  OR C X . HE  I GHT > 

DINE  NS  ION  XTV(3)jXV(3)/UV<3)»C<3/3)»B(3j3) 

C OH H ON  T ( 32B  > > IT<5 > 

ELAPSED  TIHE 

EQUIVALENCE  <T(2),TIHE) 

TIHE  STEP 

EQUIVALENCE  <T(3),DTIHE) 

ORIENTATION  HATRIX  ( BO  HB  W.R.T.  AIRCRAFT) 

EQUIVALENCE  <T(4),C) 

POSITION  OF  BOHB  C.G.  W.R.T.  AIRCRAFT  C.G. 

EQUIVALENCE  ( T < l 3 ) , X , X V > , ( T< 1 4 ) , Y > , < T<  1 5 > , 2 ) 

VELOCITY  OF  BOHB  C.G  W.R.T.  AIRCRAFT  C.G. 

EQUIVALENCE  < T< 1 6 ) , U, UV  ) , ( T<  1 7 ) , V ) , < T< 1 8 ) . W ) 

ATTITUDE  OF  AIRCRAFT  W.R  T EARTH  AXES  (EULER  ANCLES) 

EQUIVALENCE  < T < 7 6 ) , AR OL L ) . ( T < 77  ) , A P I TC H ) , < T < 78  ) , A YA W > 
BOHB  HASS  (KGH) 

EQUIVALENCE  < T ( 8 2 ) , ST HA S S ) 

GRAVITATIONAL  ACCELERATION 

EQUIVALENCE  < T( 92  ) . GRAV AC ) 
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FREE  STREAM  MACH  NUMBER 

EQUIVALENCE  ( T< 96 > , RHACHB  ) 

UNIFORM  TRAJECTORY  INDICATOR 

=-l  BEFORE  t = + l AFTER  TRANSITION  TO  BOMB-ALONE  TRAJECTORY 
EQUIVALENCE  < T< 1 14 > , FREE  ) 

FREE  STREAM  VELOCITY 

EQUIVALENCE  < T < 1 36  ) , VEL  I NF > 

GRAVITY  FORCES 

EQUIVALENCE  <T(161>,GX>,(T(162>,GY>,<T<163>,GZ> 

ORIGINAL  DIVE  ANCLE  OF  AIRCRAFT 
EQUIVALENCE  < T< 1 64  ) , D I VEB  > 

CENTRIPETAL  ACCELERATION  OF  AIRCRAFT  (IN  G'S) 

EQUIVALENCE  <T<168>,GF> 

EQUIVALENCE  <XTV(1>,XT>,<XTV<2>.YT>,<XTV<3>,ZT> 


SET  FREE  TO  +1 . 

FREE*t . 

ADJUST  X AND  U SO  THAT  THEY  NOU  REFER  TO  EARTH  AXES, 

WITH  ORIGIN  AT  GROUND  LEVEL  BELOU  AIRCRAFT  C.G.  AT  TIME  OF 
RELEASE  < T IME*B  ) . 

IF  ( ABS< GF>  GT . I . E - 6 ) CO  TO  1 
SV= VEL INF*T IME 
VGT«B 
GO  TO  2 

1 VG=2  *VELINF/( GF*GRAVAC > 

VGT  = TI ME/VG 
SV=VG*VEL INF*S IN<  VGT  ) 


2 DD=D IVEB- VGT 

D IVEaDD- VGT 

CALL  TRANS  ( X T V , AR 0 LL , A P I T CH , AY AH , - 1 . X V > 


o o o o o n o o o ui  ui  o o 
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X=XT*SV*COS(DD  > 

y *y  t 

Z*ZT  + HEI GHT 

CALL  TRANS  < X T V , AR 0 LL , A P I T CH . A V AU , - 1 , U V ) 

U = X T *V  EL  1 NF  *C  0 S(  01 VE  ) 

V * Y T 

M=ZT-VELINF*SIN(DIVE) 

C 

C REPLACE  (C)  WITH  (C).(CA),  SO  THAT  BOHB 
C ORIENTATION  IS  NOU  REFERRED  TO  EARTH  AXES. 

C 

C NOTE;  C CA  = CCA'  .C'  )' 

C CAN  TREAT  C'  AS  3 COLUMN  VECTORS  AND  USE 
SUBROUTINE  TRANS  (DON'T  USE  CTRANS) 

D031  1=1,3 
D03 1 J = 1 , 3 
B<  I , J >=C<  J , I ) 

D 03  J = 1 , 3 

CALL  TRANS  ( X T V , AR 0 LL . A P I T CH , A Y AU , - 1 , B C 1 , J ) > 

D03  1=1,3 
C( J , I )=XTV(  I > 

SET  UP  THE  NEU  DERIVATIVES 
CALL  DAUX 

CLOSE  OFF  THE  3 PROBE  CALIBRATION  FILES 

CALL  CLOSE  (3) 

CALL  CLOSE  (12) 

CALL  CLOSE  (14) 

DEFINE  THE  STORE  COEFFICIENTS  FILE 

CALL  DEFINE  ( 1 3 . 29 , 1 6 64 , ST  OR C X , I V 4 , B , B , B ) 

C 

C SET  THE  GRAVITY  FORCES  TO  EARTH  AXES  (CONSTANT) 

C 

GX=  B 
GY  = B 

GZ=-STHASS*GRAVAC 

C 

MRI TE( 9, 45) 

45  FORMATC ///l X' **BOMB-ALQNE  TRAJECTORY  B EG  I NS  ******** ** 

1EARTH  AXES  FROM  THIS  POINT************************** 
RETURN 
END 


' //  ) 
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APPENDIX  IV 

LISTINGS  OF  ARCTAN,  TRANS,  CTRANS  AND  PARAB 


C ARCTAN  DATE!  15-6-76 

C 

FUNCTION  ARCTAN<Y,X> 

FINOS  ANGLE  WHOSE  TANGENT  IS  Y/X,  IN  RANGE  B TO  TWO  PI. 

CALLED  8 Y NAIN. 

COMMON  T ( 32B  > < IT<5) 

EQUIVALENCE  < T < 1 65  > , P I > , < T ( 1 66  ) , T MOP  I ) 

C 

ARCTAN-ATAN< Y/ABS<  X ) ) 

IF  (X.LT.B.)  ARCTAN=PI-ARCTAN 

IF  < ARCTAN. LT .8.  > ARCTAN*ARCTAN  + TWOP  I 

RETURN 

END 


rv  <r>  — o n o o o o o oooooooooo  o o 
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TRANS  DATE:  4-8-76 

SUBROUTINE  TRAMS  (XT, ROLL, PITCH, YAW,  INV , X ) 

CALLED  BY  M A I N , C AL C , CRE F , F LOW 

* EULER  ANGLES  ( R 0 LL , P I TC H , Y AW  ) REFER  TO  NEW  AXES  W.R.T.  OLD  AXES. 

* ( XT,YT,ZT>=MATRIX  (EULER  ANGLES  )*( X , Y , Z > IF  INV=+1 

* (XT,YT.ZT>=  INVERSE  MATRIX  (EULER  ANGLES >*( X, Y , Z > IF  IHV»-1 

* MATRIX  ELEMENT  AIJ=CQSINE  OF  ANGLE  BETWEEN  NEW  I AXIS  AND  OLD  J AXIS. 

* I , J = 1 , 2 , 3=  X , Y , Z . 

DO  NOT  ALLOW  XT.YT.ZT  TO  BE  THE  SAME  LOCATIONS  IN  CALLING  ROUTINE 
AS  X , Y , Z , OR  IT  WILL  NOT  WORK 

DIMENSION  C(3,3>,A(3,3>,X(3>,XT(3> 

COMMON  T ( 32B>,  I T ( 5 > 

EQUIVALENCE  (T(4),C) 

EQUIVALENCE  ( A ( l , 1 ) , A 1 1 ) , ( A(  l , 2 ) . A 1 2 ) . ( A(  1 , 3 > , A 1 3 ) 

EQUIVALENCE  ( A(  2 , 1 > , A21  > , ( A(  2 , 2 ),  A 22  >,  ( A(  2,  3 ),  A2  3 ) 

EQUIVALENCE  ( A ( 3 , 1 ) , A 3 1 > , ( A( 3 . 2 ) , A 32  ) , ( A<  3 , 3 > , A3 3 ) 


SR=S IN( ROLL  ) 

CR=C0S( ROLL > 
SP=SIN(PITCH> 

CP=C0S( P I TCH  ) 
SY=SIN(YAW> 

CY*C0S( YAW) 

A 1 1 = CP  *C  Y 

A12=SR*SP*CY+CR*SY 
A13=-CR*SP*CY+SR*SY 
A 2 1 = -C  P*  S Y 
A22=-SR*SP*SY+CR*CY 
A23=CR*SP*SY+SR*CY 
A 3 i = SP 
A32=-SR*CP 
A 33  = CR  *C  P 
GO  TO  2 

ENTRY  POINT  FOR  CTRANS 

ENTRY  CTRANS  (XT,INV,X> 
DO l 1=1,3 
DGl  J= 1 , 3 
A(  I , J )=C(  I,  J ) 


o o o o o 


PARAB  DATE:  4-8-76 

SUBROUTINE  PAR A8  CABC,X,F> 

FITS  PARABOLA  F < X > = A + B *<  X -X  A > + C* < X - XA  >*  * 2 TO  DATA. 

DIMENSION  ABC<  3 > * X<  3 ) < F( 3 ) 

C 

DX2  = X(  2 >-X<  l > 

DX3  = X(  3 > - X<  1 ) 

ABCC I )=F<  1 ) 

F XA  * < F ( 3 )-F< i ) )/DX 3 

ABC( 3)=<  < F<  2)-F( 1 > >/DX2-FXA>/(DX2-DX3) 

ABC( 2>*FXA-ABC<3  )*DX3 

RETURN 

END 


ooooooooo  r>o  oooooooooo 
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LISTINGS  OF  FLOW,  CALC,  TRAV,  TRWAIT,  ADRD,  PROBE  AND  READIN 


FLOW 


DATE:  5-8-76 


SUBROUTINE  FLOU 

COLLS  CALC, TRU AIT, ADRD. PROBE, TRANS ,CTRANS 
CALLED  BY  MAIN. 

•♦♦•♦FINDS  FLOU  PROPERTIES  ALL  ALONG  BOMB  IN  BOMB  AXES 
*****AND  STORES  PROPERTIES  IN  AMS, ....UMS. 

DIMENSION  8HS<  1 2 ) , RHS<  i 2 ) , UMS(1  2 > , VMS(  1 2 ) , WMS(  12  ) , XHS ( 12),  PRESS!  5) 
DIMENSION  UPVEC<  3>,  UT<3  >,UA<  3 >,UVEC(  3> 

COMMON  T<  32B) , IT<5  > ,NF 

EQUIVALENCE  <T(96>,  RHACHB  >,  < T ( 9 9 > , AT T ACk> . <T(1B4>,  QDPB  ) 
EQUIVALENCE  ( T ( 8 3 ) , APTR OL ) , < T< 84 ) , APTP I T ) 

EQUIVALENCE  < T( 9 1 > , VSOUND  ) 

EQUIVALENCE  < T< 165  ),  PI  ) 

EQUIVALENCE  < T ( 1 93 ) , QMS ) , < T<  2B5 ) , RMS  ) , < T<  21 7 > , UMS  ) 

EQUIVALENCE  ( T< 2 29  > , VMS  > , < T<  24 1 ) , UMS  ) , < T<  25 3 > . XMS  > 

EQUIVALENCE  < T< 277  ) , DPI  3 , PRESS ) , < T< 278  > , DP24 >, 

I < T<2  79  ) ,PMAN  >, < T<  28B>, PPIT  >,<  T<  281  >, TPRESS  > 

EQUIVALENCE  (ITU),  NMS  ) 

EQUIVALENCE  < UPVEC ( 1 ) , UP  ) , ( UP VEC<  2 >, VP  > , ( UP VEC< 3 > , UP  > 

EQUIVALENCE  < UT , UVEC >,( UPVEC , UA  ) 

DATA  TPS TP/ 116  .3497/,  TP  ZERO/ 1676.  /.TPRAT/345 .3257/ 

•♦•♦♦STATEMENT  FUNCTION  FOR  TANGENT 
TAN<  X)=SIN<  X)/COS<  X ) 

DYNAMIC  PRESSURE,  1/2  RHO  V SQUARED,  IS  GIVEN  BY: 

QDPB-  (GAMMA/2  ) <RMACHB**2>  PSTATB 

UHERE  PS T AT B=  TPRESS/<  1 ♦<  GAMMA-  1 >/2RHACH8^*2  )*•<  GAMMA/C  GAMMA-  1 ) ) 
QFACT=QDPB/TPRESS 

QFACT  = 8. 7*RHACHB*^2/<  1 . + B . 2*RHACH8^2  >*^3  . 5 

IF(NFEQNMS)  GO  TO  1 

NS*  1 

NF=  NMS 

I =*NS 

I NC  = 1 

GO  TO  2 

NS=  NMS 

NF=l 


C *****  FOR  EACH  POINT  ALONG  BOMB 
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C *****  N tx  i POINT  IS  l + I NC , UNLESS  POINT  IS  LAST  POIN- 
2 NEXT  =1 +1 NC 


C *****  I F POINT  IS  LAST  POINT,  NEXT  POINT  IS  1,  WHICH  WILL  BE 
C*****CLOSE  TO  FIRST  POINT  AFTER  NEXT  INTEGRATION  STEP 
C *****  I S IT  THE  FIRST  POINT? 

IF( I NE. NS)  GO  TO  3B 

C*****CALCULATE  1ST  POINT  IN  BOMB  AXES  (METRES)  IN  TRAVERSE  AXES 
C *****  I N INCHES,  AND  START  PROBE  MOVING. 

CALL  CALC  (NS) 

C*****UAI T 'TIL  PROBE  REACHED  THIS  POINT,  AND  SETTLED 
C*****DOWN  SUFFICIENTLY  TO  TAKE  PRESSURE  MEASUREMENTS 
3 0 CALL  TRWAIT 

C*****WHEN  READY,  TAKE  PRESSURE  MEASUREMENTS 

CALL  ADRD  < PM A N , PP I T , DP t 3 . DP 2 4 , TP R ES S ) 

TPRESSM  T PR ESS -T PZERO  )/TPSTP 
D 08 4 KM, 5 

84  PRESS( K)=PRESS(K  )*TPRAT 

QDPB=TPRESS*QFACT 
C 

I F<  I .EQ.NF)  GO  TO  67 

C*****CALCULATE  NEXT  POINT  IN  TRAVERSE  AXES  (AS  BEFORE), 

C AND  START  PROBE  MOVING. 

CALL  CALC  (NEXT) 

C*****CALCULATE  FLOW  PROPERTIES  FROM  PRESSURE  MEASUREMENTS 
67  CALL  PROBE  ( D 0 WN P , S I D EP , RM S(  I ) , QM S ( I ) ) 

C*****FIND  ACTUAL  FLOW  VELOCITY  COMPONENTS  IN  PROBE  AXES  (M/SEC) 

VEL  = VSOUND*RMS<  I )*SQRT( ( 1 . +. 2*RMACHB**2  >/( 1 . 2*RMS( I >**2>  > 
UP=-VEL/SQRT( l +TAN(D0WNP)**2+TAN(SIDEP)**2> 

VP=UP*TAN<  S IDEP  ) 

WP=UP*TAN( DOWNP  ) 

C *****  F IND  VELOCITY  IN  TUNNEL  AXES  FROM  PROBE  AXES 
CALL  TRANS  < U T , A P T R OL , A P TP  I T , B . , - 1 , UPVEC) 

C*****FIND  VELOCITY  IN  PLANE  AXES  FROM  TUNNEL  AXES 
CALL  TRANS  ( U A , P I , A T T AC K , 0 . , + 1 , UT) 

C*****F!ND  VELOCITY  IN  BOMB  AXES  FROM  PLANE  AXES 
CALL  CTRAHS  (UVEC,+l,UA> 

UMS( I > = UVEC< 1 ) 

V MS ( I )=UVEC( 2 ) 

W MS ( I > = U VEC( 3 > 

C*****CONTI NUE  FOR  ALL  MEASUREMENT  STATIONS 
I F(  I . EG  NF  > RETURN 
I * I ♦ I NC 
CO  TO  2 
END 


oooooo  o o 
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CALC  DATE:  4-8-76 

SUBROUTINE  CALC  ( J ) 

••♦•♦CONVERTS  POINT  (XB,YB,ZB)  IN  BONB  AXES  (METRES)  TO 
••••POINT  (XT,YT,ZT)  IN  TRAVERSE  AXES  (INCHES) 

CALLS  CTRANS, TRANS 
CALLED  BY  FLOW 

DINE  NS  ION  XT(3),XHS(12),XA(3),X1(3) 

COHHON  T ( 32B  ) , I T ( 5 ) 

EQUIVALENCE  (T<i3>,X>,(T<14>,Y>,(T<15>,2> 
EQUIVALENCE  ( T( 8 9 ) , SCALE  >,  ( T( 99  >, ATTACK  > 
EQUIVALENCE  ( T ( 1 65  ) , PI).  (T(167>,  CNI > 

EQUIVALENCE  (T(253  ),XMS> 

EQUIVALENCE  ( XT, XI > 

EQUIVALENCE  ( X A(  1 ) , XA  1 ) , < XA(  2 ),  YA  ) , ( XA(  3 ) , Z A ) 

C 

XI<  1 )=X«S(  J ) 

X I ( 2 )=B 
Xl(  3 )=  B 

C*****TRAHSFORH  FROM  BONB  METRES  TO  PLANE  METRES. 

CALL  CTRANS  ( XA, -1 , XI  ) 

C*****THEN  CONVERT  TO  INCHES  AND  SCALE  DOWN  TO  MODEL  SIZE. 
R*CM I/SCALE 
XAla(XAl+X)*R 
YA=( YA+Y  >*R 
ZAa( ZA+Z  )*R 

C*****THEN  TRANSFORM  TO  TRAVERSE  AXES 

CALL  TRANS  (XT,  P I , AT T A C K , B . , - 1 , XA) 

CALL  TRAV  ( XT  ) 

RETURN 

END 


/TRVCTS  DATE:  8-3-77 

; 

/ROUTINE  TO  MOVE  THE  TRAVERSE  RIG  TO  A 
/POSITION  < X , Y , Z ) IN  THE  MINIMUM  TIME 
/CALLING  SEQUENCE: 

/ FORTRAN:  CALL  TRAV  (X) 

MACRO:  J MS*  TRAV 

J MP  . *2 

4BBBBB+X 

/ 

/ASSUMES  'REQUIRED  COORDS.'  ARE  SET  PRIOR 

/TO  ENTRY. 

/ 

TRSXa?051 01 
TRSY=7051 41 
TRSZ=7B5121 
T RSR*7031 61 


/ 

GLOBL  TRAV.TRWAIT 

/INTERNAL 

GLOBL  AG,  AK , AX , DA 

/EXTERNAL 

TRAV 

XX 

J MS  * DA 

J MP  + 2 

X 1 

4B0000+REQXP 

I N IT 

CAL  /SET  UP  SKIP 
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16 

/CHAIN  FOR  AN 

TRSX 

/INTERRUPT  FROM  THE 

TRINT 

/'X  COMPLETE'  FLAG 

CAL 

/.  . . 

16 

/ . . . 

TRSY 

/'Y  COMPLETE'  FLAG 

TRINT 

/.  . . 

CAL 

/ . . . 

16 

/.  . 

TRSZ 

/'Z  COMPLETE'  FLAG 

TRINT 

/ . 

CAL 

/NOTE:  THE  TRAVERSE 

16 

/SKIP  ON  ROLL 

TRSR 

/MUST  BE 

TRINT 

/INITIALIZED!  !! 

LAC  (JMP  NSU 

/DON'T  DO  THE 

DAC  INIT 

/SET-UP  AGAIN. 

NSU  DZM  XF  LG 

/CLEAR  THE  THREE 

DZM  YF LG 

/'COORDIHATE  MOVING' 

DZM  ZF LG 

/INDICATORS. 

LAC  TRAV 

/SET  UP  THE 

DAC  PCSVEi 

7B5164 

/ENABLE  TRAVERSE  FLAGS. 

LAN  -1 

/SET  THE  'TRAVERSE 

DAC  TRIP 

/IN  PROGRESS'  INDICATOR 

LAW  -3 

/SET  UP  TO  ACQUIRE 

DAC  TRKNT* 

/3  COORDINATE  VALUES. 

PXA 

/SAVE  THE  CONTENTS 

DAC  XRSVEt 

/OF  INDEX  REGISTER# 

C LX 

/AND  CLEAR  IT  . 

LAC  REQXP 

DAC  RQPTR 

T V 1 J MS  * AG 

/GET  A COORDIHATE 

ROPTR  XX 

/VALUE  X,Y/  OR  Z. 

ISZ  RQPTR 

ISZ  RQPTR 

J MS  * . AK 

/MULTIPLY  VALUE  BY 

DSA  C1BBB 

/I BBB  < DECIMAL)  . 

J MS  * AX 

/FIX  IT,  AND 

DAC  XREQ.X 

/SAVE  BINARY  VALUE 

J MS  B I NBC D 

/CONVERT  TO  BCD 

DAC  RBCD,X 

/EQUIVALENT,  AND  SAVE 

AXR  1 

/BUMP  INDEX  REGISTER. 

ISZ  TRKNT 

/ALL  3 COORDS  YET? 

JMP  TV  1 

/NO  . 

JMP  TR1 

/YES. 

i 


REOXP 
XREQ 
R B CD 
C 1 BBB 
XFLG 
YFLG 
Z F LG 


B 

BLOCK  3 
BLOCK  3 

12;  372BBB 

B /'COORDINATE  MOVING'  INDICATORS; 

fl  /NON-ZERO  IF  THE  CORRESPONDING 

B /AXIS  IS  CURRENTLY  IN  MOTION. 


.'INTERRUPT  HANDLER 

t 

TRINT  DAC  ACSV2I 
LAC*  ( B 
DAC  PCSV2* 


/SAVE  ACCR  CONTENTS. 
/SAVE  INTERRUPT 
/RETURN  ADDRESS 
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DZM  RL FL  G • 
TRSX 

ISZ  RL  FL  G 
TRSY 

ISZ  RLFLG 
TRSZ 

ISZ  RLFLG 
LAC  RLFLG 
AAC  -3 
SZA 

JMP  TRPROC 
705104 
LAC  ACSV2 
ION 

JMP*  PCSV2 


TRPROC 

LAC 

PCSV2 

DAC 

PC  S V E • 

LAC 

ACSV2 

DAC 

ACSVE# 

PXA 

DAC 

XRSVE 

TR1 

LAW 

-3 

DAC 

TRKNT 

CLX 

T R 2 

XCT 

RDTR  , X 

J MS 

TRBCD 

TCA 

TAD 

XREQ  . X 

DAC 

DEL,  X 

A XR 

1 

ISZ 

TRKNT 

JMP 

TR  2 

TRSX 

SKP 

DZM 

XF  LG 

TRSY 

SKP 

DZM 

YF  LG 

TRSZ 

SKP 

DZM 

ZFLG 

7051B4 

LAW 

-3 

DAC 

KNTT  # 

CLX 

T A 1 

LAC 

DEL,  X 

SZA 

J MP 

TA2 

AXR 

1 

ISZ 

KNTT 

JMP 

TA1 

JMP 

TRCOM 

T A 2 

LAW 

-3 

DAC 

KNTT 

CLX 

T A 3 

J MS 

TR 

AXR 

1 

I SZ 

KNTT 

JMP 

TA3 
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/TEST  WHETHER  ANY 
/OF  THE  X , Y / OR  Z 
/TRAVERSE  COHPLETE 
/'FLAGS  ARE  UP 


/ . . . 

/"RLFLG'  MILL  BE 
/3  IF  HONE  ARE  UP 

/ . . . 

/'HOT  3 LEGITIMATE  INTERRUPT. 
/'NONE  ARE  UP 
/MUST  HAVE  BEEN 
/FROM  ROLL  AXIS. 

/THEREFORE  IGNORE  IT. 

/SET  THE  RETURN 
/POINTER 

/SET  UP  TO  RESTORE 
/THE  ACCUMULATOR. 

/SAVE  INDEX  REGISTER 

/CONTENTS 

/SET  UP  TO  READ 

/3  COORDINATE  VALUES  . 

/CLEAR  INDEX  REGISTER. 

/READ  BCD  COORD  VALUE. 
/CONVERT  IT  TO  BINARY. 
/SUBTRACT  'ACTUAL' 

/FROM  ' REQUIRED' . 

/SAVE  THE  DIFFERENCE. 

/DO  THIS  FOR  ALL 
/3  COORDS  X < Y < 1 Z 
/ . . . 

/CLEAR  THE 

/"COORDINATE  MOVING' 
/INDICATORS  OF 
/THOSE  COORDS  FOR 
/WHICH  'TRAVERSE 
/COMPLETE'  FLAGS 
/HAVE  BEEN 
/POSTED 
/ . . . 

/CLEAR  TRAVERSE  FLAGS. 


/TEST  THE  THREE 
/'ACTUAL'  COORDINATES 
/ X , Y , AND  Z FOR 
/EQUALITY  WITH 
/THEIR  'REQUIRED' 

/VALUES . 

/TRAVERSE  IS  COMPLETE. 
/SET  UP  TO  EXAMINE 
/3  COORDINATES. 

/CLEAR  INDEX  REGISTER. 
/RELOAD  TRAVERSE  REGISTER 
/< IF  NECESSARY).  AND 
/BEGIN  TRAVERSING  IF 
/APPROPRIATE 
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TRRTN 

LAC 

XRSVE 

/RESTORE  THE 

PAX 

/INDEX  REGISTER, 

LAC 

ACSVE 

/THE  ACCUMULATOR, 

ION 

/TURN  ON  INTERRUPT, 

> 

J HP* 

PCS  VE 

/AND  RETURN. 

TRCOH 

DZH 

TRIP 

/CLEAR  THE  'TRAVERSE 

7 85 1 44 

/DISABLE  TRAVERSE  FLAGS. 

JHS 

TIHE 

/SET  ' TIH1 ' TO  THE  TIHE 

AAC 

6 

/AT  WHICH  TRAVERSE  COMPLETE, 

TCA 

/PLUS  6 SECONDS.  THIS  ALLOWS 

DAC 

TIH1# 

/THE  SYSTEM  TO  SETTLE. 

$ 

J HP 

TRRTN 

/RETURN . 

/ROUTINE 

TO 

TEST  WHETHER 

THE  TRAVERSE  RIG 

/INTERFACE  REGISTER  CAN  BE  FILLED  WITH 
/THE  'REQUIRED'  VALUE  OF  A COORDINATE 
/AND/OR  THE  CORRESPONDING  COORDINATE 
/CAN  8E  SET  TRAVERSING. 

/THE  ROUTINE  EXECUTES  WHICHEVER  OF  THESE 
/OPERATIONS  IT  FINDS  TO  BE  LEGAL. 


t 

TR 

XX 

LAC 

DEL,  X 

/RETURN  IF  THE 

SNA 

/'ACTUAL'  VALUE  EQUALS 

J HP  * 

TR 

/THE  'REQUIRED'  VALUE. 

DZH 

DREGR 

DZH 

TROK 

LAC 

(LAC  XFLG 

/SET  THE  SWITCHES 

DAC 

TST  1 

/IN  SUBROUTINE  'TEST' 

LAC 

(LAC  XREfi 

/SO  THAT  THEY 

DAC 

TST2 

/WILL  EXAMINE  THE 

LAC 

(TAD  DEL 

/'X'  COORDINATE 

DAC 

TST3 

/ . . . 

LAW 

-3 

/SET  UP  TO  EXAMINE 

DAC 

TRKNT 

/3  COORDINATES. 

TRBt 

JHS 

TEST 

ISZ 

TST  1 

/BUMP  THE  SWITCHES 

IS  Z 

TST2 

/IN  'TEST'  TO 

ISZ 

TST3 

/THE  NEXT  COORD . 

I SZ 

TRKNT 

/3  COORDS  YET? 

J HP 

TR  B l 

/NO  . 

LAC 

DREGR 

/OK.  TO  RELOAD 

SZA 

/TRAVERSE  REGISTER 

J HP 

TRB2 

/NO  . 

LAC 

RBCD  . X 

/RELOAD  REGISTER  IF 

JHS 

LDREG 

/'DREGR'  IS  ZERO 

TRB2 

LAC 

TROK 

/OK  TO  TRAVERSE 

SZA 

/THIS  COORDINATE? 

J HP* 

TR 

/NO 

ISZ 

XFLG,  X 

/YES  SET  INDICATOR 

XCT 

ITR,  X 

/AND  INITIATE  TRAVERSE 

J HP  * 

TR 

/RETURN . 

.’ROUTINE  TO  CARRV  OUT  THREE  TESTS  AND 
/SET  LOCATIONS  ' DREGR ' AND  'TROK' 

/ACCORDINGLY 

/UPON  RETURN  FROM  ' TEST' , 

< A ) THE  TRAVERSE  INTERFACE  REGISTER  CAN  BE 
RELOADED  I F / 6 ONLY  IF,  DREG K =8  . 

j 

u 
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( 8 THE  COORDINATE  XREQ.X  CAN  BE  SET 
TRAVERSING  IF,  t,  ONLY  IF.  TROK*B. 


TEST 

XX 

TST1 

LAC  XF  LG 

/RETURN  IF  THE 

SNA 

/COORDINATE  IS 

JNP*  TEST 

/NOT  HOVING. 

TST2 

LAC  XREQ 

/RETURN  IF  'REQUIRED' 

T C A 

/VALUES  OF  THE 

TAD  XREQ , X 

/HOVING  COORD 

DAC  TSTSV# 

/AND  THAT  OF  XREQ.X 

SNA 

/ARE  EQUAL 

JNP*  TEST 

/ . . . 

TST3 

TAD  DEL 

/RETURN  IF  XREQ.X 

X OR  TSTSV 

/IS  BETWEEN  THE 

SPA 

/HOVING  COORD  AND 

JNP*  TEST 

/ITS  'REQUIRED'  VALUE. 

ISZ  DREGR 

/DON'T  LOAD  REGISTER. 

LAC  XREQ , X 

/IS  THE  CURRENT 

TCA 

/VALUE  OF  THE 

TAD  BREGR 

/TRAVERSE  REGISTER 

DAC  TSTSV 

/BETWEEN  THE  'REQUIRED' 

TAD  DEL , X 

/COORD  XREQ.X 

X OR  TSTSV 

/AND  ITS  'ACTUAL' 

SNA 

/VALUE? 

IS2  TROK 

/NO.  DON'T  TRAVERSE. 

JNP*  TEST 

/YES.  RETURN. 

DREGR 

XX 

TROK 

J 

XX 

"R!P 

0 /'TRAVERSE  IN  PROGRESS'  INDICATOR 

T R UND 

0 

I TR 

705102  /INITIATE  TRAVERSE,  X AXIS. 

7B5142  /INITIATE  TRAVERSE,  Y AXIS. 

7B5122  /INITIATE  TRAVERSE,  2 AXIS. 

RDTR 

7B5B12  /READ 

POSITION,  X AXIS. 

7B5B52  /READ 

POSITION,  Y AXIS. 

705032  /READ 

POSITION,  2 AXIS. 

DEL 

0;  B; 

0 

DELAY  XX 

DAC  DLYKNTI 
IS2  DLYKNT 
JNP  -1 
JNP*  DELAY 
LDREG  XX 

DAC  REGR* 

J NS  TRBCD 
DAC  BREGRi 
LAW  -S065 
J NS  DELAY 
LAC  RE  GR 
705BB4 
LAN  -5 
J NS  DELAY 
JNP*  LDREG 


i 


/SAVE  THE  BCD. 

/CONVERT  TO  BINARY, 

/AND  SAVE  IT. 

/IB  HILL ISECOND 
/SETTLING  DELAY 
/LOAD  BCD  INTO 

/TRAVERSE  INTERFACE  REGISTER. 
/WAIT  FOR  REGISTER 
/TO  SETTLE  < 2B  NICROSECS), 
/AND  RETURN. 


.ROUTINE  TO  CONVERT  TO  BINARY 
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X X XXXX  XXXX  XXXX  XXXX 
S 1 8421  8421  8421  8421 


XX 

DAC 

BCD* 

DZH 

BIN* 

RTL 

DAC 

TEHP* 

SZL 

IS  Z 

BIN 

LAB 

-4 

DAC 

KNT* 

LAC 

BIN 

CLL 

HUL 

12 

LACQ 

DAC 

BIN 

LAC 

TEHP 

RTL 

RTL 

DAC 

TEHP 

RAL 

AND 

<17 

TAD 

BIN 

DAC 

BIN 

ISZ 

KNT 

J HP 

BC1 

LAC 

BCD 

RAL 

LAC 

BIN 

SPL 

TCA 

J HP  * 

TRBCD 

.'ROUTINE  TO  CONVERT  BINARY 
/TO  8421  BCD. 

('CALLING  SEQUENCE: 

.'WITH  THE  BINARY  IN  THE  AC  , 
/JHS*  BINBCD 
/NORMAL  RETURN 
-'BCD  IS  RETURNED  IN  THE  AC  . 


8 I MB  CD 


BB  1 


XX 

DAC  BB  IN* 
DZH  BCD* 

DZH  SIGN* 
SNA 

J HP  BB  1 
TCA 

DAC  BB  I N 
LAC  (4BBBBB 
DAC  SIGN 
LAU  -S 
DAC  BKNT* 
LAC  BB  IN 
LMQ 

SPL ! CLL ! CLA 

LAC  ( B4BBBB 

DAC  LINK* 


6 B 2 
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C LA 
0 IV 
12 

XOR  BCD 
RTR 
RTR 

XOR  LINK 
D AC  BCD 
I SZ  BKNT 
J HP  BB2 
RAL 

XOR  SIGN 
J HP  * BINBCD 

$ 

/ 

/ROUTINE  TO  READ  THE  HACRQ  TINE  FROH  THE  REAL  TINE 
/CLOCK,  AND  CONVERT  TO  SECONDS  AFTER  HIDNIGHT. 

/ 

/CALLING  SEQUENCE: 

/ JHS  TINE 

/ NORHAL  RETURN 

/TIHE  IS  RETURNED  IN  THE  AC. 

t 

TIHE  XX 

7B2512 
CHA 

JHS  TAUX 
JHS  HUL6B 
DAC  TRTIHE 
7B2552 
CHA 
DAC  HINSEC 
RCL 
SWHA 

AND  <377 
JHS  TAUX 
TAD  TRTIHE 
JHS  HUL6B 
DAC  TRTIHE 
LAC  HINSEC 
AND  <377 
JHS  TAUX 
TAD  TRTIHE 
J HP  * TIHE 

/ 

TRTIHE  B 
HINSEC  B 


/ 

TAUX  XX 

DAC  TEHPX 

RTR 

RTR 

AND  <17 
CLL 
HUL 
12 

LACQ 

DAC  TE  HP  X* 1 
LAC  TEHPX 
AND  <17 
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TAD  TENPX+1 
J HP  * TAUX 

Jt 

TEMPX  0 ; B 


MUL6B  XX 
CLL 
HUL 
74 

LACQ 

J HP  * HUL6B 

t 

/ROUTINE  TO  WAIT  AT  LEAST  6 SECONDS  FROM  THE  TINE 
/OF  COMPLETION  OF  THE  PREVIOUS  TRAVERSE. 

/ THE  CURRENT  TINE  IS  READ  FROH  THE  REAL  TINE  CLOCK 
/ AND  IS  COMPARED  WITH  THE  CONTENTS  OF  LOCATION 
/ ' T I M 1 ' « WHICH  WAS  SET  AT  'TRAVERSE  COMPLETE'. 

/THE  PROGRAM  IS  TRAPPED  UNTIL  BOTH  THE  'TRAVERSE  IN  PROGRESS' 
.'INDICATOR  (TRIP)  IS  CLEAR  AND  THE  CURRENT  TINE  IS  GREATER 
/THAN  THE  CONTENTS  OF  TIN!. 

/CALLING  SEQUENCE: 

/ J MS  * TRWAIT 

NORMAL  RETURN 

TRWAIT  XX 

TW1  LAC  TRIP 

SZA 

JNP  TW1 
LAC  TI Ml 
SNA 

JMP*  TRWAIT 
T W 2 J MS  TIME 

TAD  TIN! 

SPA 

JMP  TW2 
DZM  TIN1 
JMP*  TRWAIT 

END 


,'ADRD  DATE:  21-4-Zb 

$ 

/ROUTINE  TO  READ  THE 

MANIFOLD  PRESSURE 
PITOT  PRESSURE 

DIFFERENTIAL  PRESSURES  DP13  t DP24 

/ 

FROM  THE  RAYTHEON  HULTIVERTER, 

'AND  THE  SI  TOTAL  PRESSURE  FROH  THE  DIGITIZER 

# 

.'Calling  SEQUENCE:  CALL  ADRD  <PHAN,PPIT,DP13.DP24,TPRESS> 


.’CALLED  BY  FLOW 


& T = i 
hD=6 


H 1 


hddvex 


h 2 


HODCH 
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IODEV 

DT,  AD 

GLOBL 

ADRD 

GLOBL 

DA,  AG, . AH 

/ 

GLOBL 

CPBIN 

HDRD 

XX 

J NS  * . 

DA 

JNP  +6 

PN 

a 

PP 

a 

P 1 3 

a 

P 2 4 

a 

TP 

1 

a 

ONCE 

INIT 

AD, B, ADRD 

INIT 

DT , B, ADRD 

SEEK 

DT , CONPBL 

READ 

DT , B, ANPVEX, 

WAIT 

DT 

CLOSE 

DT 

LAC  VEX1P3 
T AD  (TAD  ADBUF 
DAC  ADDVEX 
LAC  (JNP  A 1 
DAC  ONCE 


AU 

AN 


AH 


READ  AD, B, ADBUF, 16 
WAIT  AD 
LAC  ADBUF+22 
TCA 
XX 

J NS  * 

J NS  * 

C5 

J MS  * 

V EX  5 
PXA 

DAC  XRSVE* 

PLA 

DAC  LSVE# 

cl:: 

LAC  (4 
PAL 

LAC  (ANPGNS 
DAC  CAINP 
LAC  ( PCAL 
DAC  PC  ON  S 
LAC  PM , X 
DAC  PPTR 
LAC  ANPBL.X 
TAD  (TAD  ADBUF 
DAC  AD  DC  H 
LAC  AH  P2  E R , X 
TCA 
XX 

J NS  * AU 
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J MS*  . AL 
GAINP  XX 

I S2  GAINP 
IS Z CA IMP 
J MS*  AK 
VEX  5 

J MS*  .AH 
TEMP 

J NS  POLY 
PCONS  XX 

TEMP 

PPTR  XX 

LAC  PCOHS 
AAC  6 
DAC  PCONS 
AXS  1 
JMP  A2 

/ 

/READ  TOTAL  PRESSURE 

» 

t 

7B3412  /READ  S1TP. 

J MS  * CPS  I N 
J MS  * .Ail 
LAC  TP 
DAC  TPTR 
<J  MS  * AH 
TPTR  XX 

LAC  <5 
PAL 

2ARG  DZH  PH , X 

AXS  1 
JMP  ZARG 
LAC  XRSVE 
PAX 

LAC  LSVE 
PAL 

JMP*  ADRD 

/ 

/DATA  AREA: 

Jt 

ADBUF  BLOCK  23 

/ 

TEMP  0 ; B 

V E X5  a;  B 

CS  3;  24BBBB  /S . 

COMPBL  SIXBT  'C0MPBLDA8' 

/ 

/DATA  INFORMATION  BLOCK: 

J 

AMPVEX  174BBBJ  B 
DATSCE  . BLOCK  2B 

AMPBL  BLOCK  12 

NAMP  B 

VEX1P3  B 

CHKNT  B 

GP2PR  B 

AMB  BLOCK  2B 

VEXBL  BLOCK  2B 

SYSREQ  BLOCK  S 

MACREQ  B 


o o o o 


RtIPGNS 

BLOCK 

24 

ampzer 

BLOCK 

12 

PCAL 

. BLOCK 

74 

THERMS 

. BLOCK 

6 

THERMC 

BLOCK 

6 

BALK 

BLOCK 

1 10 

CPRES 

BLOCK 

4 

/CALLING  SEQUENCE: 

JUS  POLY 

PTR  TO  CAL  CONSTS 

PTR  TO  VOLTAGE 

/ PTR  TO  PRESSURE  RETURN. 

/ NORMAL  RETURN 


/ 

POLY  XX 

LAC*  POLY 
ISZ  POLY 
D AC  PL  2 
AAC  2 
0 AC  PL  3 
AAC  2 
DAC  PL  5 
LAC*  POLY 
ISZ  POLY 
DAC  PL  1 
0 AC  PL  4 
LAC*  POLY 
ISZ  POLY 
DAC  PL  6 


J MS  * 

AG 

PL  1 

XX 

J MS  * 

. AK 

/V 

P L 2 

XX 

J MS  * 

A I 

/A  * V 

P L 3 

XX 

J MS  * 

. AK 

/A* V+B 

P L 4 

XX 

J MS  * 

. AI 

A A *V  +B  ) * V 

P L 5 

XX 

J MS  * 

AH 

/A . V2+B  . V + C 

P L 6 

XX 

J 

JHP* 

POLY 

/ 

END 

ADRD 

C PROBE 


DATE:  31-10-77 


USES  DIRECT  ACCESS  I/O. 

SUBROUTINE  PROBE  ( DONNP . SI DEP . RMACH. Q ) 

C 

C CALLS  READIN. 

C CALLED  BY  FLOW  . 

C 

DIMENSION  DOWNS( 2)<  SIDESC2  ). ANGLE<  2) 

COMMON  T ( 320).  IT<5  > 
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EQUIVALENCE  ( T ( 1 34  ) , D TR  ) , ( T(  1 35  ) , R TD  ) 

EQUIVALENCE  < T(  277  ) , DPI  3 ), < TC  278  ) , DP24  ) 

EQUIVALENCE  < T< 2 79  ) , PHAN  ) . < T<  28B > - PP I T > , < T<  281  ), T PRESS  ) 
C 

EQUIVALENCE  ( ANGLE(  1 ) , DOWN  ),  ( ANCLE( 2 ), S IDE  ) 

C 

DATA  ERROR/B. 15/.PITCHN/B. / 

C 

TAN<  A )*S I N(  A ) / COS( A > 

C 

C*****FIND  PRESSURE  RATIO 


PRATIO=PP IT/PHAN 
IF  (PRATIO.GT. 1.  > GO  TO  43 
WRI TE( 4, 42  > PPIT.PMAN 

42  F0RMAT(/1X' (PROBE)  P I TO T = ' , G 1 2 . 5 , ' < MANIFOLD*'  G12.5) 
STOP  5 

C *****  S ET  P I TCH=*L  AST  PITCH  FOUND  ( I . E . . AT  LAST  MEASUREMENT 
C*****STATI ON/ WHICH  VALUE  WILL  BE  CLOSE  TO  THIS  ONE) 

43  PITCH=PITCHN 
C*****GO  THRU  MAX  OF  5 TIMES 

NTHRU=B 
NTMAX1 =9 

C«*«**F  IND  MACH  NUMBER  AT  THIS  RATIO. PITCH 
IB  IF  ( NTHRU  LE  4 ) GO  TO  44 

WRITE( 4/ 45)  PITCH, DIFF.PRATIO.RMACH. DOWN, SIDE 
45  F0RMAT(/1X' (PROBE)  5 ITERATIONS  EXCEEDED'  / 

1 IX'PITCM. DIFF, PRATIO, RMACH, DOWN. SIDE'  / 

2 1X6C12  5) 

STOP  11 

4 4 IR=6  5* IF IX( PITCH  )+ I F I X( 2B  *( PRAT  I O-l  . ) > + l 


CALL  READINL 14. RMACH. 1.2. PRATIO. PITCH. G5.  IR  ) 
AMI =RM ACH  *1 B 
N 1 = A M 1 -3  . 

IF  ( Ml  GE  1 AND  Ml  . LE  . 1 1 > GO  TO  53 
WRITE( 4. 52)  RMACH 


52  F0RMAT(/1X' (PROBE)  MACH  NO.  * '.G12  5.'  OUTSIDE  RANGE') 
STOP  12 

C*****FIND  TOTAL  PRESSURE  AT  THIS  MACH, PITCH 

53  IR=1 1 * I F IX( PITCH  )*M1 

CALL  READIN(3. RATIO, 1,2, RMACH, PITCH, 11, IR) 

PTOT=PPIT/ RATIO 
C***»*F  IND  STATIC  PRESSURE 

PSTAT=PTOT/’((  1 +B.2*RMACH**2)**3.5) 

C****»F  IND  DYNAMIC  PRESSURE- 

Q=B5*1.4*PSTAT*RMACH**2 
C*»***F  IND  COEFFICIENTS 
DCl 3*DP1 3/Q 
DC24*DP24/Q 

C • ****F  IND  MACH  NUMBER  INTEGERS  ON  EACH  SIDE  OF  ACTUAL  MACH  NO 


C • ** **F OR  EACH  MACH, FIND  DOWN-AND  SIDE-WASH  ANGLES  (DEGREES) 

IR=(  (Ml  — 1 )*41*IFIX(  lfl.*(DC24*2.  ) ) ) *4  1 ♦ I F I X(  1 B . •(  D C 1 3*  2 > >♦  1 
1F(  IR  LE  10*16  81  ) GO  TO  47 
WRITE( 4, 46)  IR. RMACH, DC13. DC24 


44  F0RMAT(/1X' (PROBE)  RECORD  NO.  *',I6,'  EXCEEDS  11*1681'/ 

1 lX'RNACH, DC13, DC  2 4 : '.3G15.6) 

STOP  1 
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20  l R* I R+  168  1 

C****«L INEAR  INTERPOLATION  BETWEEN  MACH  NUNBERS 
F MACH* Ml 

FHACH*AH l -FHACH-3 . 

OOWHPaDOWNS( 1 >+( D0WNS(2 >-DOWNS( 1 >>*FNACH 
S IOEP*SI DES( 1 >♦( SI DES<  2 )-S IDES< 1 ) )*FNACH 
C*****C0NVERT  PROBE  ANGLES  TO  RADIANS 
DQWNPaDOWNP*DTR 
SIDEP-SIDEP*OTR 

C*»**»C ALCULATE  NEW  PITCH  IN  DEGREES 

P I TCHN  = AT  AN  (SORT  (TAN  < DOWNP  )**2*TAN  < SI DEP >**2  ) )*RTD 
C*****ABSOLUTE  DIFFERENCE  BETWEEN  NEW  AND  OLD  PITCH  • 

DIFF*ABS  < P I TCHN -P I TCH  > 

C*****SET  OLD  PITCH  EOUAL  NEW  PITCH- 
P ITCH=PI TCHN 

C***»* INCREASE  NO  OF  TINES  GONE  THRU 
NTHRUaNTHRU*l 

C * ** ** C OH PA R E WITH  ERROR  SET  ABOVE 
IF  ( DIFF . GT  ERROR  ) CO  TO  IB 
C*****RETURN  WITH  PROBE  ANGLES  IN  RADIANS 
RETURN 
END 


READIN 


DATE:  15-7-76 


SUBROUTINE  READIN  < ND AT , F. NF , IND, X , Y . IXD I H. I R ) 

C 

C CALLED  BY  CREF, PROBE. 

C 

C * **  **  T H I S SUBROUTINE  READS  TABULATED  FUNCTIONS  FROM  DISK 
C*****FROH  . DAT*2 . AND  PERFORMS  A 2-D  LINEAR  INTERPOLATION  TO  FIND  THE 
C*****VALUES  OF  THE  FUNCTIONS  AT  THE  POINT  (X.Y). 

C*****CALL  IT  WITH  THE  FOLLOWING  P ARAHETERS . . . 

C • ****F( 1 > ...  F<  NF  )*  UP  TO  IB  INTERPOLATED  FUNCTION  VALUES  (RETURNED). 
C*****NF*NUHBER  OF  FUNCTIONS  TO  BE  INTERPOLATED. 

C*****IND-NO  OF  INDEPENDENT  VARIABLES  IN  RECORD  (REAL  ♦ INTEGER) 
C*****X,Y-POINT  AT  WHICH  WANT  INTERPOLATED  VALUE  OF  FUNCTION. 

C*»***DATA  STRUCTURE  ON  DISK  ASSUMED  TO  BE  OF  THE  FORM... 

C • ** ** 1 ST  RECORD  - NX , NY  (DIMENSIONS  OF  X.Y  TABLES) 

C *****  2N0  RECORD  - ( X T ( I ) , I * 1 , NX  ) ( YT( I ), I -1 . NY  ) TABULATED  VALUES. 
C*****NEXT  NX*NY  RECORDS  - EACH  OF  THE  FORM 
C *****  I X IY  (12)  XT(  I X > YT(IY)  (ZT(IZ))  F(  1 ) ....  F(HF) 

C*****WITH  IX  VARYING  HOST  RAPI DLY  < THEN  IY  (THEN  12. IF  PRESENT). 
C*****INTERPOLATES  BY  USING  2X2  POINTS  ABOUT  ACTUAL  ONE. 

C*****I  E . . THE  4 TABULATED  POINTS  CLOSEST  TO  ACTUAL  POINT  (X.Y). 

C 

C*****X,Y  TABLES  (HAXINUH  OF  65  VALUES  AT  PRESENT) 

DIMENSION  XT( 2 . 2 ) 

C • •• **UP  TO  IB  FUNCTION  VALUES  CAN  BE  READ  IN 

C*****UP  TO  3 INTEGER  AND  REAL  NUNBERS  IX , I Y . I Z . XT( IX ) . YT< I Y > . ZT< I Z > 
DIMENSION  I T( 3 ) , RT( 3 ) 

C*****2X2*4  POINTS  FOR  EACH  OF  UP  TO  II  FUNCTIONS  WANTED  FOR  INTERP 
DIMENSION  FW( II. 2.2  ) 

C *****THE  INTERPOLATED  VALUES  (UP  TO  IB)  ARE  RETURNED  THRU  F 
DIMENSION  F ( 1 ) 

COMMON  T ( 321).  1 1 T( 5 > 

EQUIVALENCE  (T(135).RTD) 
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DO  1 BB  1*1/2 
DOIBB  J=1 ,2 

IRN=IR  + IXDIH*<  J-l >♦< 1-1  > 

READtMDAT'  IRM  )<  IT<  L )/ L*t  , IMO  )/<  RT<  L)/L“l/  IND  )» 
l < F H(  L < I , J ).  L«1  /HF  ) 

C IFCNDAT.EQ  13)HRITE(9,3B1)  ( FM<  L<  I < J )»L  = l»NF)/NOAT.  IRN  » 

C 1 < 1 T(  L > / L*  1 / IND) 

C 3B1  F0RHAT<1X'FW  ( R/ I H > ' / IX 8G 1 5 . 6 / 1 X2G 13 . 6 , ' NDAT/1RN.1T 
J K=2-I ABS< I -J > 

IBB  XT(  1/J  )=  R T<  JK  > 

C IF<NDAT.EQ  1 3 > WR  I T E < 9 , 3 B B ) ( ( XT<  I , J ) / J ■ 1 , 2 ) , I*  1 , 2 > 

C 3BB  F OR  N AT< lX'XT  < READ  I H > ' . 4G1 5 . 6 > 

C *****  C ARRY  OUT  2-D  INTERPOLATION 
C*****FOR  EACH  FUNCTION 

DX=<  X-XT<  1, 2>  XT<  2.  1 >-XT<  l , 2)  > 

DY=< Y-XT<  1, 1 > >/<  XT<  2, 2>-XT< 1 , 1 > > 

D012B  N= 1 i NF 

C ***** 2 -D  INTERPOLATION, AND  RETURN  RESULT  THRU  F 
C*****FOR  THE  POINT  FU(N.l,J> 

F 1 1 = FU < N / 1,1) 

F 12  = FU<  N / 1/2) 

F21 =FH<  N / 2/  1 ) 

F 22  *Fy < N / 2/ 2) 

12B  F<N  )=F11*( 1 ,-DX-DY  )+F21*DX+F12*DY  + <Fll+F22-F12-F21  )*DX*DY 

RETURN 
END 
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APPENDIX  VI 

LISTINGS  OF  CORECT  AND  CREF 

C CORECT  OR TE  : 26-7-76 

C 

SUBROUTINE  CORECT  ( ALPHA. U > 

C 

C CALLS  PARAB 

C CALLED  BY  MAIN 

C 

C*  THIS  ROUTINE  CONSTRUCTS  THE  STREANLINE  ABOUT  THE  REFERENCE  POINT  AND 
C*  THUS  CALCULATES  THE  CHORD  ANGLE  ALPHA  AND  CORRECTION  ANGLE  GANNA 

C*  THIS  CONSTRUCTION  IN  THE  Z-X  PLANE  IS  BASED  ON  THIS  ANALYSIS-  ((RADIAN 

C*  SLOPE  OF  STREANLINE  AT  POINT  X « DZ/DX  = ( DZ/DT >/(  DX/DT > 

C*  = U(X)/U(X>  (IN  TERNS  OF  STREAH  VELOCITY  CONPONENTS). 

C*  HENCE  THE  STREANLINE  IS  DEFINED  BY  Z( X )* I NTEGR AL  FRON  ZERO  TO  X 

C*  OF  U<  X >/U(  X ) . 

C 

DINE  NS  ION  U(12),H(l>,X(i2>,XV(3>,F(3>,ABC(3> 

CONNON  T ( 32H  > , IKS  ) 

EQUIVALENCE  (IT<2>,  NREF  ) 

EQUIVALENCE  <T(217>,  U>,  (T(  253>,  X> 

EQUIVALENCE  < ABC(  l ) , A >,  ( ABC(  2 >,  B > . ( ABC(  3 > , C > 

EQUIVALENCE  < X V< 1 > , XA  ) , ( XV < 2 > , XB ) , ( X V( 3 > , XC  ) 

•FIND  THREE  POINTS  ABOUT  THE  REFERENCE  POINT* 

I XA*NREF- 1 

IF  ( NREF . EQ . 1 ) IXA=1 
IX*1XA 

•DEFINE  POINTS  XA,XB,XC  ABOUT  NEASURENENT  STATION  XB* 

INTEGRAND  F = U/U 
DOl  L*  1 , 3 

XV( L > = X< IX)-X(  IXA) 

F<  L )««( IX  >/U<  IX) 
l IX*IX+1 

C • F I T PARABOLA  F=A+B*X+C*X**2  TO  THE  INTEGRAND* 

CALL  PARAB  (ABC, XV, F> 

C * I NT  EG  R ATE  IT  FRON  ZERO  TO  EACH  POINT  XA,XB,XC  TO  OBTAIN 
C *STREAHLINE  POINTS  ZA,ZB,ZC,  THUS  - Z*A*X+B*X**2/2*C*X**3/3  * 

ZBX  B *<  C/3  . * XB  + B/2 . )*XB+A 
Z CX  C *( C/3  . * XC  + B/2 . ) *X C* A 

C *STRAI GHT-L I NE  CHORD  DEFINED  BY  Z*( ZC/XC  )*X . * 

C *D  IFFERENCE  BETBEEN  CHORD  AND  STREANLINE  AT  REFERENCE  POINT* 

D I F F*( ZBXB-ZCXC  )*XB 

C****» ALPHA  IS  CHORD  ANGLE,  TUOGAN  IS  CORRECTION  ANGLE. 

ALPHA= AT  AN( ZCXC  ) 

TUOGAN  *2  *D IFF /XC 

C 

C CALCULATE  EFFECTIVE  INCIDENCE 

C 

ALPHA* ALPHA+TWOGAN 

RETURN 

END 
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C CREF 
C 
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DATE!  1-7-77 


SUBROUTINE  CREF 

C 

C CALLS  TRANS , READ  IN . 

C CALLED  BY  MAIN.  [SYSTEM 

C*  CONFUTES  FORCE  AND  NOMENT  COEFFICIENTS  DUE  TO  WIND,  IN  TOTAL  AXIS 
C*  BY  INTERPOLATION  INTO  DATA  TABLES 
C*  PITCHT  <8  2 4 .38)  ROLLT  (-45  -37.5. .45) 

C*****RMACHT  ( .4,  .5, .6, .7.  .8,  .85,.?,  .95) 

C *****  THE  COEFFICIENTS  ARE  ON  DISK 
C*****EACH  RECORD  WITHIN  A FILE  CONTAINING 

C*****I.J,K»PITCHiROLL.HACH,CX,CY,CZ,CL,CH,CN,CLP,CHQ,CNP,CYP . 

C 

DIMENSION  C(1B),CINT(2,18),PTVEC<3) 

C * **  **  (CINT(M,I  ) CONTAINS  COEFF  I AT  TWO  MACH  NOS  H> 

DIMENSION  UMS(  12  ).  VMS(  12  >,  WMS(  12).  RNS(  12) 

DIMENSION  RNACHT(8> 

COMMON  T ( 32B), IT(5  ) 

EQUIVALENCE  < T ( 9 4 ) , TH ET 0 T ) , < T ( 95 ) , PH  I T 0 T ) 

EQUIVALENCE  (T(98),DMAX) 

EQUIVALENCE  < T ( 1 29  ) , P ) , < T<  1 38  ),  Q ) , < T<  1 3 1 ) , R ) 

EQUIVALENCE  < T(  l 34  ) , DTR  ) , < T< 1 35  ), RTD  > 

EQUIVALENCE  (T(165),PI> 

EQUIVALENCE  (T(285),RMS) 

EQUIVALENCE  < T < 2 1 7 ) , UMS  ) , ( T<  229  > , V MS  ), ( T< 24 1 > , WHS  ) 

EQUIVALENCE  ( T< 265  ) , CX  ) , ( T< 286  ) . CY  ), < T<  267 ) , CZ ) 

EQUIVALENCE  ( T < 2 68  ) , CL  ) , < T ( 2 6 9 ) , C M ) , < T < 27B ) , CN ) 

EQUIVALENCE  < IT(2),NREF) 

EQUIVALENCE  < C ( 1 ) , T CX  ) , ( C<  2 ) , TC Y ) , ( C < 3 ) , T CZ  ) 

EQUIVALENCE  ( C( 4 >, TCL  ), ( C<  5 > , TCM ) , < C< 6 ) , TCN  ) 

EQUIVALENCE  ( C ( 7 ) , T CL P ) , < C< 8 > , TCM Q , T CN R ), < C ( 9 ) , T C NP  ) 

EQUIVALENCE  <C(1B),TCYP) 

EQUIVALENCE  < PTVEC(  1 ) , PT  ), < PT VEC<  2 ). QT  ) , < PTVEC( 3 ) , RT ) 

C*****MACH  NUMBERS  AT  WHICH  COEFFS  ARE  TABULATED 
DATA  RMACHT/.4,  .5,  .6,  .7,  .8,. 85.  .9,  .95/ 

C 

THDEG=THETOT*RTD 

C 

IF< THDEG . GE  8 . AND . THDEG  LE. 38.  ) GO  TO  3 
WRITE<  4,232  > THDEG 

232  FORNAK/1X'  THETOT  ='  F12.2,'  OUTSIDE  RANGE') 

STOP  3 

3 PHIDEG=PHITOT*RTD 

C***»*GET  ROLL  (PHI)  FROM  (8  TO  368)  INTO  RANGE  OF  TABLES  (-45  TO  *45) 

C 

I PH  I =(  PH  I DEG+45  . )/9B  . 

PHIDEG  = PHIDEC-98.*FL0AT(  IPHI  ) 

C*****FlND  MACH  NUMBER  INTEGERS  ON  EITHER  SIDE  OF  ACTUAL  ONE 
RMSREF  = RMS( NREF  > 

D 0 1 HA CH  = 1 , 8 

IF(RMSREF  GE  RHACHT(HACH). AND. RMSREF.LT. RMACHT(HACH  + 1 ))GO  TO  2 
1 CONTINUE 
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y R I TE( 4/21)  RdSREF 

21  F0RMAT</1X'<CREF ) MACH  HO*  ',G12.5*'  OUTSIDE  RANGE') 

C * E XPRE S S EULER  RATES  IN  TOTAL  AXIS  SYSTEM* 

2 CALL  TRANS  < P T * P H I T OT  * 0 . , 0 . , - 1 , P > 

C *****  F OR  EACH  MACH  NUMBER  (MACH  AND  HACH+1) 

IR=<  MACH-1  >*208+ IF  I X<  (PHIDECM5.  )/7. 5 >•  16*1  F IX<  TKDEG/2  >♦! 
DO600  N=l*2 

C WRITE< 9, 4291 > IR * MACH * PH  IDEG * THDEC 

C 4291  FORMATC 1 X' IR, MACH, PHIDEC. THDEG' * 16, 16- 2G15  6 ) 


C*****INTERPOLATE  TO  GET  10  COEFFS  AT  THIS  PITCH*  ROLL. MACH. 

CALL  READIN( 1 3 < C * 10*3*  THDEG*  PHIDEG. IS*  IR> 

C *****  S TORE  10  COEFFS  AT  THIS  MACH  NO 
DO  10  1=1*10 
10  C INT<  N * I )=C(  I ) 

C*****D0  FOR  BOTH  MACH  NUMBERS*  EACH  SIDE  OF  ACTUAL  ONE 
b 0 0 I R=  I R+  20  8 


C*****NOH  INTERPOLATE  BETWEEN  MACH  NUMBERS  TO  GET  FINAL  COEFFS 
RMACH1 =RMACHT< MACH  ) 

RHACH2«RMACHT ( MACH+1  ) 

DO  500  1=1*10 

SL0PE  = <CINT<2.  I )-CINT< 1 * I ) )/( RHACH2-RMACH1  ) 

5B0  C< I >»CINT< 1 * I >♦<  RMSREF-RMACH1  )* SLOPE 

C * A DD  COMPONENT  PARTS  TO  GIVE  RESULTANT  COEFFICIENTS* 

VELREF=SQRT(UMS<  NREF >**2+VMS<  NREF )**2*WMS< NREF  )**2) 

C T CX  = TCX 

TCY=TCY+TCYP*PT*DMAX/2. /VELREF 
C TC2=TC2 

TCL=TCL+TCLP*PT*DMAX/2. /VELREF 
TCM=TCM+TCNQ*QT*DHAX/2. /VELREF 
TCN  = TCN+<  TCNP*PT-TCNR*RT  >*DM AX/2 . /VELREF 
C*****TRANSFORM  TO  BOMB  AXES  FROM  TOTAL 

CALL  TRANS  C CX , P H I T OT * 0 . * 0 . , + 1 . TC X > 

CALL  TRANS  ( C L * P H I T OT  * 0 . * 0 . * + 1 . TC  L ) 

RETURN 

END 
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APPENDIX  VII 
LISTING  OF  CINC 

C CINC  DATE : 26-8-76 

C 

SUBROUTINE  CINC 
CALLED  BY  HAIN 

CALCULATES  INCREHENTS  IN  COEFFICIENTS  DUE  TO  FLOW  VARIATION  ALONG  BONB 
IN  BONB  AXES. 


DINE  NS  ION  ASTOR<  12),D(12),Q(  1 2 > , R<  12),U(12),V(  12).W(12),X<  12) 
DINE  NS  ION  F YL  < 12  ),  F ZL  < i 2 ),  FNL  ( 12)  , FNL(  1 2 ) , F YZ(  12 , 2 ) , DF<  2 , 2 ) 

DINE  NS  I ON  D ADX( 1 2 ) , DVDX<  12),DWDX<  12),0DX<  12.3) 

D INENS ION  ABCC3) 

CONNON  U 32B  > . I T < 5 > 

EOU  I VALENCE  < T<  9 1 ) . VS  >,  < T<  97  ),AHAX  >,  <T<  98),  DHAX  ) , ( T(  1 04  ) , ODPfl  ) 
EQUIVALENCE  ( T < 9 6 ) , RN ACHB ) , < T< 1 36 ) , VEL I NF ) 

EQUIVALENCE  < T < l 69  ) , A ST OR > , < T< 1 81 ) , D >, ( T< 1 9 3 > . Q ) 

EQUIVALENCE  ( T < 2 1 7 > , U >,  < T<  229  ) . V ) , < T<  24  1 > , W ) 

EQUIVALENCE  <T<253  ),X) 

EQUIVALENCE  < T<  2 72  > , DCY  ) , ( T<  273  ),  DC2  ),  < T<  275  ),  DCH  >,  < T<  276),  DCN) 

EQUIVALENCE  ( I T(  1 ) , NNS  ) , ( I T<  2 ),  NREF  ) , ( I T<  3 > , NB  ),  < I T<  4 ) . NF  ) 

EQUIVALENCE  ( DADX. DDX< 1 , 1 ) ) , < DVDX , DDX< 1 , 2 ) ) , ( DWDX , DDX( 1,3)) 
EQUIVALENCE  < ABC< 1 > , A > , < AB C< 2 ), B > , < ABC< 3 > , C ) 

EQUIVALENCE  ( F Y2  ( 1 , 1 ) , F YL  , FNL  ),  < F Y2<  1 . 2 ),  F2L  , FNL  > 


DATA  ETACDC/B . 864/ 


SUNF<DX ) = (<  C*DX/3  . + B/2.  )*DX+A  )*DX 
SCSFCANG  >=SIN( 2. *ANG)*COS(  ANG/2  . ) 

NFN2=NF-2 


FIND  DA/DX,  DV/DX  AND  DW/DX. 

NETHOD: 

IF  F = A ♦ B( X-XA  > ♦ C<  X-XA  >**2 , 

THEN  DF/DX  = B + 2C<X-XA). 

C IF  THE  DERIVATIVE  IS  REQUIRED  AT  POINT  J.  THEN  TNE 

C ROUTINE  FITS  PARABOLAS  SUCCESSIVELY  THROUGH  THE  SETS  OF  POINTS 
C (J-2,J-I,J).  <J-1,J,J*1>  AND  < J , J* 1 , J +2  ) , TO  OBTAIN 
C 3 ESTINATES  OF  DF/DX  BY  THE  FORHULA  ABOVE. 

C PROGRAH  THEN  AVERAGES  THE  THREE  ESTINATES.  AT  THE  END 

C POINTS.  THE  TREATNENT  IS  SLIGHTLY  DIFFERENT. 

C 

D04  L* 1 , 3 
DOI  1=1, NNS 
1 DDX<  I , L > = B 

HBN2-NHS-2 
D02  I * 1 , N BN  2 
X A=  X ( I ) 

GO  TO  (9,  IB, 1 t ),  L 

9 CALL  PAR AB  < AB C , X<  I ) . AS T OR < I ) ) 
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GO  TO  12 

IB  CALL  PAR  AB  (ABC.  X<  I >.  V<  I > ) 

GO  TO  12 

11  CALL  PARAB  < ABC  , X(  I > » W<  I ) > 

12  DDX(  I , L > = DDX<  l , L ) + B 
11*1+1 

D0X< I1.L>*DDX< I1,L)+B+2.*C*(X<I1>-XA> 

12=1+2 

2 DDX< 12 , L )=DDX(  1 2 , L ) +B +2 . *C*< X( I 2 >-XA  ) 

IF  ( NHS . EQ . 3 ) GO  TO  4 
DDX  ( 2, L >=DDX<  2,L >/2  . 

DDX<NNS-l,L)“DDX(NHS-l,L>/2. 

IF  ( HNS. EQ. 4>  GO  TO  4 
D03  1=3, HBH2 
DDX(  I » L ) = DDX<  I , L )/3  . 

CONTINUE 

CONPUTE  VALUES  AT  REFERENCE  POINT. 

UR=U( NREF  ) 

VR=V<  NREF  ) 

W R=W( NREF > 

QR=Q( NREF  ) 

VISCOUS  FORCE/DC  I)  = 1/2  RHO  ETA  CDC  V SQRT<V»*2+W**2  ) 
WHERE  RHO  = 20/ ( U**2 + V** 2 +W**2 > . 

R1  = QR*ETAC0C*S0RT(  VR**2  + WR**2  )/< UR**2  + VR**2  + WR**2  ) 

V YR  = R 1 *V  R 
VZR=R1 *WR 

SLENDER  BODY  F OR CE/( D A/DX > = Q S1H<2  ANG  > C0S(AHG/2) 
WHERE  ANG  = ARCTAN<W/U>  FOR  FY 
AND  = ARCTAN<  V/U  > FOR  FZ  . 

ALPHA= AT AN( WR/UR  ) 

BETA«ATAN< VR/UR) 

OAL PH=QR*SCSF < ALPHA  ) 

QBE  T *QR*SCSF< BETA) 

C 

C AT  EACH  HEASURENENT  STATION,  CALCULATE  THE 

C CURVED  FLOW  INCRENENTS. 

C 

D05  1=1, NHS 
QI  = Q< I ) 

U1  = U<  I ) 

VI  = V< I ) 

W I = W( I ) 

RHO *2  *QI/<  UI**2+VI**2+WI**2  ) 

C 

C SLENDER  BODY  INCRENENTS. 

C 

ALPHA= AT  AN( WI /UI  > 

BET  A=ATAN( V I/U  I ) 

DFYSB=<QI*SCSF<BETA  )-QBET  )*DADX( I ) 

DFZSB=<OI*SCSF< ALPHA >-QALPH)*DADX<  I > 

C 

C BUOYANCY  INCRENENTS. 

C 

C B'JOYANCY  FORCE  = - VELINF  RHO  A DV/DX  FOR  FY,  AND 
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- VELINF  RHO  A DV/DX  FOR  FZ. 
THUS.  IN  UNIFORM  FLOW.  BUOYANCY  FORCE  IS  ZERO. 


Rl=VELINF*RHO*ASTOR< I ) 

DFYB»-R1*DVDX<  I ) 

DFZB>-Rl*DUOX( I ) 

VISCOUS  INCREMENTS. 

Rl-B  5*RH0*ETACDC*SQRT<  VI**2  + UI**2  ) 
DFYV«<R1*VI-VYR>*D<  I) 

DFZV«<R!*UI-VZR>*0<  I) 

F YL ( I > = DFYSB+DFYB*DFYV 
3 FZL< I)*DFZSB+DFZB+DFZV 

C NRITE<  9.  123  )FYL.  FZL 

C 123  F0RNAT(/1X' FYL '/1X1BG13 . 5/1X2C13. 5// 

C 1 lX'FZL ' /1X1BC13. 5/1X2C13 .5//) 

C 

C INTEGRATE  U.R.T.  X FROM  X(NB>  TO  X(NF>.  TO  GET  TOTAL 
C FORCE  AND  MOMENT  INCREMENTS. 

C 

D08  1=1,2 
DO?  J= 1 , 2 
FF=B  . 

K=NB 

6 CALL  PARAB< ABC.XCK >,FYZ< K, J > ) 

FF=FF  + SUMF<  X<  K+2  >-X<K  >> 

K =K  ♦ 2 

IF<K.LE.HFM2)  GO  TO  6 

IF  < K. EO . NF  ) GO  TO  7 

CALL  PARAB( ABC ,X(NFM2  ), FYZ<NFM2,J  > > 

FF=FF  + SUNF< X< NF  )-X< K> > 

? DF<  I , J )=FF 

D 08  L» 1 , NMS 
FNL<  L )*X< L )*FYL<  L ) 

8 FML( L)  = -X(L  >*  F ZL  < L > 

OAM=QDPB* AMAX 
OAHD=QAM*DMAX 
DCY  = DF( 1 , 1 )/Q AH 
D CZ  » DF  < 1 , 2)/0AM 
DCM* DF( 2 , 2 >/0 AMD 
DCN  = DF( 2 , 1 >/Q AMD 
RETURN 
END 
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APPENDIX  VIII 
LISTINGS  OF  INTM  AND  DAUX 


C I NTH  DATE:  26-7-76 

C 

SUBROUTINE  INTN 
CALLED  BY  MAIN. 

* INTEGRATES  18  SIMULTANEOUS  FIRST  ORDER  DIFFERENTIAL  EQUATIONS. 

* REQUIRES  D I HENS  I ON  OF  T TO  BE  4»N  + 3 C*75> 

*** 

*****N*N0  OF  EQUATIONS  TO  BE  INTEGRATED 

COMMON  T< 32B  > / IT<5  ) 

EQUIVALENCE  C TC 2 > . T IME > / ( TC 3 > , DT > 

ENTRY  POINT  FOR  INTEGRATING  ONE  TIME  STEP. 

* CALLED  BY  MAIN. 

FOURTH  ORDER  RUNGE  KUTTA  METHOD 

T < 1 )-TIME 
DOl  J = 4 . 21 

1 TCJ+36>=TCJ> 

DO  6 J =3 , 6 
A*C  8-J  )/2 
B*J/2 

12  DO  66  1*4/21 

KN= I +36 
ICO*  I +54 
KP* I +1 8 

IF  < J . EQ . 6)  GO  TO  2 
T<  K 0 )*T(KO>  + DT«TCKP  >*B 
T C I ) = TC  KN > + DT*T( KP  >/A 
GO  TO  66 

2 TCI  ) = TCKN)+CDT*TCKP>+T<K0>>/6. 

TCKO  )* 0 . 

66  CONTINUE 

TIME-TC 1 )+DT/A 
CALL  DAUX 
6 CONTINUE 

RETURN 
END 
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C DAUX  DATE : 7-7-77 

C 

SUBROUTINE  DAUX 
C 

C CALLED  BY  INTM. 

C 

C*  ROUTINE  SPECIFIES  THE  EQUATIONS  OF  NOTION  IN  FIRST  ORDER  FORN  SO 

C*  THEY  CAN  BE  INTEGRATED  ONE  TINE  STEP  BY  INTH 

C*  FORCE  EQUATIONS  IN  AIRCRAFT  AXES. 

C*  HOHENT  EQUATIONS  IN  BOMB  AXES. 

C 

DIHENS  ION  C(3.3).X(3).U(3).H<3).ABC< 3, 6 >. TINES(3  ) 

D I HENS  ION  DC<3.3).DX(3>.DU(3).DH(3) 

DINENSION  FX(3).FL(3).A<3).P(3> 

DIHENS  ION  ARQT( 3 . 3 ) . BROT( 3 . 3 ) < ACC( 3 ) 

COHHON  T < 32B  > . IT<5  > 

EQUIVALENCE  <T<2>.TIHE> 

EQUIVALENCE  < T<  4 >.  C >.  < T<  13  >.  X ).  < T<  It  >.  U >.  < T<  19  >,  H > 

EQUIVALENCE  < T(  2 2 ) . DC  >.  < T<  31  > , DX  > . < T(  3 4 >»  DU  ) . < T<  37  > , DH  ) 
EQUIVALENCE  < T ( 7 9 ) . A ) , < T< 82 >. STHASS > 

EQUIVALENCE  < T< 11 1 ) , PA  ) . ( T< 1 1 2 ) , QA  >. < T< 1 1 3 > , RA  ) 

EQUIVALENCE  < T < l 1 4 > » F RE  E > 

EQUIVALENCE  <T<117>.ACC> 

EQUIVALENCE  < T ( 1 2B > . F X ) , < T < 1 2 3 ) . F L > 

EQUIVALENCE  < T < 1 29  > , P ) , < T< 1 3B  ) , Q ) , < T< 1 3 1 ) , R > 

EQUIVALENCE  ( T< 1 37  ) , AROT  ) . < T< 1 4 6 ) , BROT  ) 

EQUIVALENCE  < T < 3 BB  ) , T I HE S ) , < T ( 3B3  ) , ABC > 

C 

DELT  *T IHE-T IHES<  1 > 


C*  RENORHAL ISE  DIRECTION  COSINES 
DOl  1*1.3 
SUH*0 

UPDATE  ANGULAR  VELOCITIES  FRON  ANGULAR  HOHENTUH  AT  SANE  TINE 

P<  I ) =H<  I )/A<  I ) 

D 02  J * 1 . 3 

2 SUH  = SUH+C( I .J  )**2 

SUH  = SQRT<  SUH  ) 

DOl  J = 1 . 3 

1 C ( I , J )*C(  I. J ) /SUH 

C*  UPDATE  BOHB  AND  AIRCRAFT  ROTATION  HATRICES 


IF  < FREE . GT  .8  > GO  TO  7 
ARO T< l .2  )=R A 
AROT< l .3  > = -QA 
AR0TC2. 1 )*-RA 
AR0T(2>3  >*PA 
A RO  T ( 3 . 1 ) *Q  A 
AROT<  3 » 2 ) *- PA 
7 B RO  T<  1 . 2 > *R 

BROT( 1 .3  ) *- Q 
BRO  T( 2 . 1 > = - R 
BROT <2.3  >=P< 1 ) 

BRO  T( 3 . 1 )*Q 
BROT ( 3 . 2 )*-P<  1 ) 

C*  EQUATION  OF  BOHB  ORIENTATION  U.R.T.  AIRCRAFT. 

C*  DC*BROT*C-C*AROT 

C 
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EQUATION  OF  BOHB  C.G  NOTION  REFERRED  TO  AIRCRAFT  AXES. 
XA«CA(  XE-XE* ) 

DXA*CA( DXE-DXE*  > + DCA( XE-XE*  > 

BUT  DCA= AROT . CA 
THEREFORE 

DXA»CA( DXE-DXE* >+AROT ,CA( XE-XE*  ) 

LET  U=CA( DXE-DXE* ) 

THEN 

DXA=U+ AROT  XA 
AND 

DU  = C AC D2XE -D 2XE* >+AROT  CA( DXE-DXE*  ) 

I E . 

DU“CA . D2XE-CA  D2XE*+AR0T. U 

WHERE 

CAD2XE*F0RCE S/HASS 
CA.D2XE*=AIRCRAFT  ACCELERATIONS 

«(  CFG  SIN( ATTACK), 0. -CFG . COS( ATT ACK  > ) 

DX“U+  AROT*X 

FORCE  EQUATIONS  OF  NOTI ON , REFERRED  TO  AIRCRAFT  AXES. 
DU=FX/STHASS-ACC+AROT*U 

NONENT  EQUATIONS  OF  HOT  ION , REFERRED  TO  80HB  AXES. 

DH*FL  + BROT*H 

DOS  1 1 = 1,3 
DOS  1 J = 1 , 3 
SUN  * B . 

DOS  K= 1 , 3 

SUH»SUH+BRQT(  I ,K  )*C(K, J ) 

IF  (FREE  LT  B.)  SUH*SUH-C<  I , K >*AROT<  K,  J ) 

CONTINUE 
DC<  I , J >=  S UH 
D 00  1=1,3 
DX<  I )=U<  I ) 

DU< I )■<<  ABC(3, I )*DELT+ABC(  2, I >>*DELT+ABC( 1, I >>/STHASS 
IF  ( FREE . LT  . 0 . ) DU<  I ) = DU( I )- ACC< I ) 

13=1+3 

DH< I )=< ABC<  3,  I3)*DELT+ABC<  2, 13)  )*DELT  + ABC< 1 , 13  ) 

D 06  J=  1 , 3 

IF  (FREE . GT  0 . ) GO  TO  6 
DX(  I >*DX< I)+AROT(I,J)*X(J> 

DU< I )=DU(  I > + AROT( I , J )*U( J ) 

DN( I )=DH( I )+BROT( I , j )*H( J ) 

URI TE( 9,44)  D X , DU , DH 

F0RHAT(/1X' DX , DU , OH  (DAUX)  ',9C12.5//> 

HRI TE( 9,46)  S TNASS , DELT . H 
F0RHAT(1X‘STHASS,DELT,H  1 ,5612.5/ > 

URI  TE(  9,  47  ) ( ( ABC(  I , J ),  1=1 ,3  ),  J = 1 . 6) 

FORH AT( 1 X ' ABC  ' , 3C 1 5 . 6/( 5X3C 1 5 . 6 ) > 

RETURN 

END 
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Figure  1.  Definition  of  real  world  coordinate  systems 
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Figure  2.  Definition  of  wind  tunnel  coordinate  systems 
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Figure  3.  The  camber  of  a curved  stream 
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U r - | V | COS  0T 
V =-  |V)  SIN  0 7 SIN  -y 
ZB  W = - | V | SIN  0f  COSy>T 


XB>XT 

Figure  4.  Definition  of  Wind  Vector  (V)  Relative  to 
Total  (T)  and  Bomb  (B)  axes 


TAN(DOWNP):  Wp  / Upl 

TAN  (SIDEP)  = Vp/llp  I BY  DEFINITION 

THEREFORE: 


Figure  5.  Definition  of  Wind  Vector  (V)  Relative  to  Probe  axes 
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